Anopheles bionomics in a malaria endemic area of southern Thailand

Background Ivermectin mass drug administration (MDA) could accelerate malaria elimination in the Greater Mekong Subregion. This study was performed to characterize the bionomics of Anopheles in Surat Thani province, Thailand. Methods Mosquitoes were collected via human landing collections between February and October 2019. Anopheles mosquitoes were morphologically identified to species. Primary Anopheles malaria vectors were dissected to assess parity status, and a subset were evaluated for molecular identification and Plasmodium detection. Results A total of 17,348 mosquitoes were collected during the study period; of these, 5777 were Anopheles mosquitoes. Morphological studies identified 15 Anopheles species, of which the most abundant were Anopheles minimus (s.l.) (87.16%, n = 5035), An. dirus s.l. (7.05%, n = 407) and An. barbirostris s.l. (2.86%, n = 165). Molecular identification confirmed that of the An. minimus s.l. mosquitoes collected, 99.80% were An. minimus (s.s.) (n = 484) and 0.2% were An. aconitus (n = 1), of the An. dirus (s.l.) collected, 100% were An. baimaii (n = 348), and of the An. maculatus (s.l.) collected, 93.62% were An. maculatus (s.s.) (n = 44) and 6.38% were An. sawadwongporni (n = 3). No Anopheles mosquito tested was Plasmodium positive (0/879). An average of 11.46 Anopheles were captured per collector per night. There were differences between species in hour of collection (Kruskal–Wallis H-test: χ2 =  80.89, P < 0.0001, n = 5666), with more An. barbirostris (s.l.) and An. maculatus (s.l.) caught earlier compared to An. minimus (s.l.) (P = 0.0001 and P < 0.0001, respectively) and An. dirus (s.l.) (P = 0.0082 and P < 0.001, respectively). The proportion of parous An. minimus (s.l.) captured by hour increased throughout the night (Wald Chi-square: χ2 = 17.31, P = 0.000, odds ratio = 1.0535, 95% confidence interval 1.0279–1.0796, n = 3400). Overall, An. minimus (s.l.) parity was 67.68% (2375/3509) with an intra-cluster correlation of 0.0378. A power calculation determined that an An. minimus (s.l.) parity reduction treatment effect size = 34%, with four clusters per treatment arm and a minimum of 300 mosquitoes dissected per cluster, at an α = 0.05, will provide 82% power to detect a significant difference following ivermectin MDA. Conclusions The study area in Surat Thani province is an ideal location to evaluate the impact of ivermectin MDA on An. minimus parity. Graphical abstract

Between 2010 and 2019, the number of malaria cases in the GMS fell by 90%, and the target of the countries within the GMS is to eliminate malaria by 2030 [1]. Malaria transmission in the GMS is complex, with various Anopheles species in the Dirus complex, Minimus complex and the Maculatus group having been recognized as the primary malaria vectors in the region. From the Minimus complex, Anopheles minimus (s.s.) is the primary malaria vector and can be found across Thailand and the GMS [2][3][4][5][6]. There are two members of the An. dirus complex that are primary malaria vectors; of these, An. dirus (s.s.) occurs east of the Thai-Myanmar border, and An. baimaii occurs west of the Thai-Myanmar border [2,6,7]. Additionally, members of the An. maculatus group are considered to be primary vectors that contribute to malaria transmission, including An. maculatus (s.s.) [2,4,5] and to a lesser extent An. sawadwongporni [2,4,6], with the latter identified as a malaria vector in southern Thailand [2]. Members within a species complex differ in their behavioral characteristics, which in turn drives Plasmodium transmission dynamics; therefore, an accurate identification of mosquitoes is essential to design and evaluate vector control methods in the GMS.
The diversity of vector species, insecticide resistance and increasing antimalarial drug resistance are some of the greatest challenges for malaria elimination in the GMS [8]. Moreover, the outdoor-feeding, outdoor-resting and early-evening feeding behaviors of GMS malaria vectors [9][10][11] also limit the effectiveness of indoor residual spraying (IRS) and long-lasting insecticide-treated nets (LLINs). Therefore, novel vector control measures which target these outdoor-feeding vectors could accelerate malaria elimination in the GMS. The mass drug administration (MDA) of ivermectin to humans has been suggested as a possible malaria parasite transmission control tool as it directly kills Anopheles that feed on treated people, regardless of blood-feeding time or location. Evidence from West Africa has shown that a single ivermectin MDA can kill wild An. gambiae (s.l.) [12,13], shift the population age structure [13], reduce the Plasmodium falciparum sporozoite rate [13,14] and that repeated ivermectin MDAs reduce clinical falciparum malaria episodes [15]. In the GMS, at human-relevant concentrations, ivermectin is lethal to An. dirus, An. minimus, An. sawadwongporni and An. campestris, and it inhibits the development (i.e. sporogony) of Plasmodium vivax in An. dirus and An. minimus [16]. A clinical trial demonstrated that An. dirus fed blood from persons treated with ivermectin (400 µg/kg) within 6 days post ivermectin administration had a 50% reduced survival rate; the reduction in An. minimus survival increased to 90% when they were fed blood from treated persons up to 10 days post drug administration [17]. This evidence has inspired a largescale cluster randomized trial in Thailand to assess the effect of ivermectin MDA on entomological and parasitological parameters of malaria transmission.
Evaluating vector control interventions for malaria in the GMS is difficult as transmission occurs primarily in the forest due to Anopheles habitat preference, combined with sporadic human entry into the forest for various agricultural and economic pursuits, not all of which are legal [18,19]. However, most vector control interventions (e.g. LLINs and IRS) are applied in the village setting where transmission is less likely to occur, and thus they do not directly target the areas of active Plasmodium transmission. Thailand has dramatically reduced its malaria burden and has set the goal to eliminate malaria by 2024. Indeed, from 2000 to 2019, there was a 96% reduction in number of malaria cases, from 159,120 to 5832, but the ratio of Thai to non-Thai cases has increased from 57 to 72% demonstrating that there are still active foci of transmission in Thailand [20]. This reduction in malaria complicates evaluation of Anopheles vector control interventions in Thailand due to low rates of Plasmodium transmission. To evaluate the effect of ivermectin MDA, a site in Thailand needed to be selected where persons live and work with active malaria transmission. Rubber plantations in Thailand offer an ideal location for Anopheles vector control intervention evaluation, as mature rubber plantations tend to be located in hilly areas adjacent to natural forests, both of which are ideal habitats for primary malaria vector proliferation. Rubber tappers live and work in the same environment, and tappers work throughout the night exposed to wild Anopheles; consequently, rubber tappers have higher rates of malaria than non-rubber tappers [21][22][23]. Indeed, a seven-province wide survey of case data from malaria clinics in Thailand from 2013 to 2016 found that 60.1% (3330/5541) of P. falciparum cases were identified from rubber tappers [23].
Of all the provinces in Thailand, Surat Thani province has the largest rubber plantation coverage, approximately 3829 km 2 [24], and is one of the few provinces affected by malaria that is not along an international border. Interestingly, in Surat Thani, the dominant malaria species infecting humans is P. falciparum, accounting for 77.58% (519/669) of all malaria cases from 2015 to 2019 [25]. Plasmodium falciparum is the ideal parasite to assess during vector control interventions as it is most sensitive to transmission interruption due to its non-relapsing nature. A cross-sectional molecular malaria survey conducted in Surat Thani in 2019 demonstrated that persons who stayed outdoors during the night-time were at the highest risk of malaria infection [26]. In that study, three districts in Surat Thani province, namely the Phanom, Vibhavadi and Khiri Rat Nikhom districts, were selected for evaluation based on their malaria case incidence comprising 65.13% (338/519) P. falciparum cases between 2015 and 2019 and high rubber plantation coverage. Historically, IRS with 5% deltamethrin has been performed in villages with a higher malaria case burden, but this control measure ceased in Phanom and Vibhavadi in 2015 and in Khiri Rat Nikhom in 2016. LLINs are widely distributed throughout all three districts by government and non-governmental organizations. Daytime indoor thermal fogging for Aedes aegypti control with 1% deltamethrin or 25% cypermethrin still occurs focally and sporadically in response to local dengue cases (personal communication, Surat Thani Vector-borne Diseases Control Center 11.3, Surat Thani, Thailand).
Little information regarding Anopheles vector bionomics in Surat Thani has been published. The largest study conducted in Surat Thani to date reported a total of 3778 Anopheles mosquitoes collected from February 2015 to December 2016 via human landing collection (HLC) in Phanom district [27]. In the same study, six Anopheles species were collected, including An. dirus (s.l.), An. minimus (s.l.), An. maculatus (s.l.), An. barbirostris (s.l.), An. hycranus (s.l.) and An. tessellatus. The predominant species was An. minimus (s.l.), comprising 87.19% of Anopheles collected. The highest mosquito densities were found between March and May in both years [27]. Molecular identification has verified the presence of An. maculatus (s.s.) in Phanom district [28] and An. minimus (s.s.) in Khiri Rat Nikhom district [29]. Due to active P. falciparum transmission, asymptomatic malaria observation, presence of primary malaria vectors and an ideal environment for evaluating a vector control intervention, Surat Thani province was selected as the study area to evaluate the impact of ivermectin MDA on entomological and parasitological parameters. Due to decreasing rates of Plasmodium transmission in Thailand, entomological (i.e. mosquito population age structure) and parasitological (i.e. human malaria prevalence) outcomes will be used to assess ivermectin MDAs. The unexpectedly low prevalence of malaria in the study sites [26] necessitates a strong emphasis on the former.
Baseline entomological surveillance utilizing the mosquito HLC method were performed to evaluate vector abundance, composition, landing activity and parity rates linked to molecular identification in Surat Thani in 2019. These efforts were undertaken to determine whether the study area was appropriate to evaluate the impact of ivermectin MDA on An. minimus population age structure (i.e. parity).

Description and maps of field sites Field site description
Collection sites were located in Khiri Rat Nikhom, Phanom and Vibhavadi districts, in Surat Thani province (651 km south of Bangkok). Villages with higher malaria case burden, as reported by local Vector Borne Disease Units, were selected from each district. Populations and house locations were mapped (see section Mapping process), and cluster sizes of approximately 300-500 persons were established. In total, 13 clusters were selected for entomological evaluation: five clusters in Khiri Rat Nikhom, four clusters in Phanom and four clusters in Vibhavadi. Khiri Rat Nikhom district is located in the center of Surat Thani province (9°1′48″N, 98°57′12″E), with its western part located in the forested hills of the Phuket mountain range adjacent to Khao Sok National Park and its eastern part mostly consisting of flat terrain. Phanom is located in the southwest of the province (8°51′18″N, 98°48′48″E) and is covered by mountains and forest. Its northwestern part is protected by the Khao Sok National Park and its southwestern part is protected by the Khlong Phanom National Park. Vibhavadi is a small district situated in the north-central portion of the province (9°14′20″N, 98°58′44″E) and is covered by mountain and forest. The western part of the district is protected by the Kaeng Krung National Park and Khlong Yan Wildlife Sanctuary. In Surat Thani, the dry season occurs from January to February, and the rainy season lasts from March to September; heavy monsoon rains occur October to December.

Mapping process
Latitude and longitude coordinates of all houses and HLC locations were captured using a 60CSx GPS unit (Garmin, Olathe, KS, USA). The open-source QGIS software was used to generate maps for each cluster. The terrain data were derived from a topographic map with elevation contour lines at 25 m.

Mosquito collection
Adult mosquitoes were collected from February to October 2019 for two consecutive nights per cluster per month, with the exception of Vibhavadi where collections began in June. Mosquitoes were collected by the HLC method. Mosquito collector volunteers were local Thai residents, non-pregnant and non-breastfeeding adults (age 18-62 years) who were capable of providing informed consent and capable of comprehending the HLC method. A canopy was constructed at each outdoor collection site to protect the mosquito collector volunteers from the elements. There were two collection sites per cluster each night. Collection sites were chosen based on proximity to a forest or rubber plantation, presence of a house nearby for access to water and electricity and close proximity to potential An. minimus (s.l.) larval habitat. Efforts were made to sample as many areas in each cluster that met the above criteria and were safely accessible. If a collection site yielded few or no Anopheles, then the site was switched to a new location the following night. At each collection site, two mosquito collector volunteers worked together from 18:00 h to 24:00 h and then replaced by two other collectors from 24:00 h to 06:00 h. Collections occurred for 50 min of each hour, followed by a 10-min break. Volunteers were instructed to wear double-layer clothing. Each collector exposed only their legs and captured mosquitoes as soon as they landed on them using a plastic collection tube sealed with a cotton ball. The mosquitoes were then transferred to collection cups grouped by hour of collection and separated by each volunteer. The cups were kept in a Styrofoam box and covered with moist towel to keep the mosquitoes humid and alive during transport back to the field station.

Mosquito morphological identification and parity dissections
At the field station, mosquitoes were transferred to plastic knockdown chambers and anesthetized with triethylamine (Flynap ® ; Carolina Biological Supply Co., Burlington, NC, USA) for 2 min, following which they were identified morphologically under a dissecting stereomicroscope (Stemi 305; Carl Zeiss AG, Oberkochen, Germany) using a standard key of adult Anopheles of Thailand [30]. All primary Anopheles species, such as the An. dirus complex and An. maculatus group, and a subset of 20 An. minimus complex per collector pair per site from 18:00 h to 24:00 h and from 24:00 h to 06:00 h were dissected to remove their ovaries. Once mosquitoes were anesthetized, they were placed in plastic Petri dishes and care was taken to keep the mosquitoes in an ice-chilled cooler with a damp towel to maintain humidity. Thus, almost all mosquitoes were alive at the point of ovary dissection. Ovaries were dissected with minuten pins in a drop of bottled Crystal water, a brand of marketed water commonly found throughout Thailand. The last two abdominal segments were gently pulled apart from the abdomen, and the ovaries were separated from the remaining internal organs and transferred to an individual well on a 12-well slide. The slide was then allowed to air dry, and care was taken to observe each ovary before total evaporation of water to assess parity status. Determination of parity status was based on the presence of coiled (nulliparous) and uncoiled (parous) tracheole skeins viewed at ×10 and ×40 magnification with a compound microscope (model B-190TB; Optika Srl, Ponteranica, Italy) and images of the ovaries were taken for reference. On occasion, if the primary reviewer had doubts on the parity status, then a second reviewer was consulted. The mosquito sample was then bisected between the thorax and the abdomen and stored in labeled 2-ml centrifuge tubes with silica gel desiccant. Processed mosquitoes were shipped back to the Armed Forces Research Institute of Medical Sciences in Bangkok for molecular species identification and Plasmodium infection evaluation of the thorax.

Molecular methods for Anopheles and Plasmodium identification DNA extraction method
The DNA extraction method involved adding 700 µl of phosphate buffer saline (pH 7.4) and 4.5-mm steel beads (Copperhead; Crosman Corporation, Bloomfield, NY, USA) to a 2-ml tube containing an individual mosquito thorax or abdomen, followed by homogenization in the TissueLyser II (Qiagen, Hilden, Germany) at 22 Hz for 2 min. The mosquito suspension was then centrifuged at 12,000 rpm for 5 min, and a 250-µl aliquot of the supernatant was used for DNA extraction according to the QIAsymphony ® DNA Minikit and Tissue LC 200 DSP protocol in the fully automated QIAsymphony ® SP system (Qiagen). The DNA was eluted in 50 μl and stored at − 20 °C until further use. DNA/RNA-free distilled water was included in the extraction process as a negative extraction control.

Molecular methods for Anopheles sibling species identification
To identify anopheline sibling species, multiplex allelespecific PCR (AS-PCR) assays were used to examine the internal transcribed spacer 2 (ITS2) region of DNA [33]. Previously published protocols [31][32][33][34] were used, with the following modifications. The amplification was carried out using total volumes of 25 μl, with the final optimized reaction conditions as follows: (i) 1× GoldStar Best Master mix (GoldStar DNA Polymerase, dNTP, PCR stabilizer and enhancer); (ii) three specific primer cocktails, each containing four or five different primer pairs to discriminate between species with 400 nM for the primers specific for the Funestus group and Maculatus group and 500 nM for the primers specific for the Dirus complex; (iii) 4% dimethyl sulfoxide (DMSO) was included only for Dirus complex reactions; (iv) a universal forward primer, located in the conserved region of the 5.8S gene, and species-specific reverse primers in the ITS2 spacer region were employed to amplify a portion of the mosquito ITS2 region; and (v) 1 μl of genomic DNA was used as template. Positive mosquito controls [i.e. An. minimus (s.s.), An. dirus (s.s.), An. sawadwongporni]) were obtained from mosquito colonies maintained at the Armed Forces Research Institute of Medical Science Department of Entomology. Negative DNA/RNA-free distilled water controls were included.
Amplifications were performed in a T100 DNA thermal cycler (Bio-Rad, Hercules, CA, USA) under the following PCR conditions. For the Maculatus group, the cycling program started with an initial denaturation at 95 °C, 10 min; then denaturation at 94 °C/1 min, primer annealing at 55 °C/30 s, extension at 72 °C/30 s for 35 cycles; and a final extension at 72 °C for 10 min. The amplification conditions for the Dirus complex were the same as those for the Maculatus group, except that the annealing time was 15 s instead of 30 s. For the Funestus group, which contains the Minimus complex, the cycling program started with an initial denaturation at 95 °C, 10 min; then amplification at 94 °C/30 s, 45 °C/30 s, 72 °C/40 s for 35 cycles; and a final extension at 72 °C for 10 min. The amplified PCR products were subjected to DNA fragment analysis using the QIAxcel Advanced System (Qiagen) according to the manufacturer's instructions. Briefly, 10 μl of PCR product was analyzed with the QIAxcel DNA Fast Analysis Cartridge (Qiagen), using the DM190 method and QX 15-bp/1-kb alignment markers. Fragment sizes were calculated using the BioCalculator (Qiagen).

DNA sequencing and data analysis
To confirm the results of multiplex AS-PCR assays, a representative of each anopheline group was selected for confirmation with the ITS2 rDNA gene in the DNA sequencing assays. Amplification of the ITS2 rDNA gene from mosquito DNA extracts was conducted using universal primer ITS2A (5′-TGT GAA CTG CAG GAC ACA T-3′) and ITS2B (5′-TAT GCT TAA ATT CAG GGG GT-3′) [35,36]. Reactions were performed in a T100 DNA thermal cycler (Bio-Rad). PCR reaction mixtures (25 μl) consisted of 2 µl of mosquito DNA extract, 0.1 U of AmpliTaq ® Gold DNA Polymerase (Life Technologies, Thermo Fisher Scientific, Carlsbad, CA, USA), 1× Gold buffer, 0.2 mM of dNTP, 2 mM of MgCl 2 and 0.2 µM of each primer. The PCR cycling program consisted of an initial denaturation step at 94 °C, 10 min; then 94 °C/30 s, 60 °C/1 min, 72 °C/1 min for 37 cycles; and a final extension step at 72 °C for 5 min. The size of the PCR product was determined using the QIAxcel Advanced System as described above. The PCR product was cleaned by the ExoSAP-IT ™ PCR Products Clean-up kit: 2 µl of ExoSAP-IT ™ was added directly to 5 µl of PCR reaction product, incubated at 37 °C for 15 min and then at 80 °C for 15 min. The ITS2 rDNA gene PCR product was sequenced using the Bigdye ® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA, USA) according to Applied Biosystems' protocol with forward and reverse universal ITS2A and ITS2B primers and run on a SeqStudio Genetic Analyzer (Life Technologies). Forward and reverse nucleotide sequence data for each sample were assembled using the Sequencher 5.1 software package (Gene Code Corp., Ann Arbor, MI, USA). The Anopheles species were identified by phylogenetic analysis. Briefly, the MUSCLE algorithm was used for sequence alignment in Molecular Evolutionary Genetics Analysis (MEGA) 6.0 software [37]. Maximum likelihood trees were constructed with the best fit model of nucleotide substitution with bootstrapping (1000 replicates) using MEGA 6.0 software as described by Tamura et al. [37].

Real-time PCR for detection of Plasmodium
Real-time PCR for Plasmodium detection from Anopheles thoraxes was performed using a 7500 Fast Real-time PCR System (Life Technologies). The primers and probes were modified from a previously published protocol of Kimura et al. [38]  cycles; the cut-off values were set automatically at every run. All positive samples were confirmed for Plasmodium species by nested PCR.

Nested PCR Analysis for Plasmodium species differentiation
Nested PCR was performed using two amplification processes as described by Kimura et al. [38] to identify four Plasmodium species: P. falciparum, P. vivax, P. ovale and P. malariae. The nested PCR was performed in a T100 DNA thermal cycler (Bio-Rad). For the first PCR reaction, the reaction mixture (20 µl Nuclease-free water was used as a negative control. The amplified fragment was analyzed using the Qiaxcel Advanced System with an expected size of around 140-160 bp. For the second PCR, the product of the first PCR reaction was diluted (1:50) with nuclease-free water and used as a template in the second PCR reaction. The reaction mixture of the second PCR reaction (20 µl) and the cycling program were the same as those of the first PCR reaction, except that the reverse specific primers were used instead of the P2R primer. The expected product (about 110 bp) of each Plasmodium species was determined using the Qiaxel Advanced System as described above.

Statistical analysis
Mosquito diversity per cluster was estimated using the Shannon-Wiener index [39] and the Simpson index [40]. The number of Anopheles caught per night was calculated by assessing the total number of Anopheles captured at one HLC site divided by two for the pair of collectors used at each location; this number was then used to calculate the mean number of Anopheles. The median mosquito catching time was calculated, and the comparison of the distribution of mosquito catching times among groups were assessed using the non-parametric Kruskal-Wallis H-test; pairwise comparisons for significantly different groupings were performed using the post-hoc Dunn's multiple comparisons test to account for the multiple comparisons using Prism version 7.2 (GraphPad Software, San Diego, CA, USA). Parity by cluster and hour was assessed by a logistic regression model, and the clustering of outcomes was accounted for by using the robust standard errors. The sample size and power calculations for a cluster randomized trial were performed using the clustersampsi command in STATA version 14 (StataCorp, College Station, TX, USA).

Topographical maps
A terrain map of Surat Thani province with the clusters evaluated from Phanom (PN), Khiri Rat Nikhom (KR) and Vibhavadi (VB) districts is shown in Fig. 1. Figure 2 is a zoomed-in version of Fig. 1, with the terrain maps depicting the house and mosquito collection locations for each of the 13 clusters evaluated in this study.
The Shannon-Wiener index and Simpson index were calculated per cluster for all species collected and for only the Anopheles species collected. For all species captured, no clear difference in mosquito diversity was found across the three districts, with the four clusters (PN-01, PN-04, KR-09, VB-10) of highest diversity occurring in each of the three districts. However, for Anopheles species, the four Phanom clusters had the highest diversity ( Table 1).
The average number of Anopheles specimens collected per volunteer per night is shown in Table 2. Overall, nightly HLCs did not vary by district for all Anopheles, with an average of 11.46 Anopheles captured per collector per night. There were no differences in mean number of An. minimus (s.l.) captured per collector per night across the three districts. However, more An. dirus (s.l.) were captured per person per night from In Phanom, 12 Anopheles species were morphologically identified from four clusters. Anopheles minimus (s.l.) was the most abundant Anopheles species captured across all four clusters (Fig. 5a-d), and the second and third most abundant Anopheles species by cluster were An. dirus (s.l.) and An. barbirostris (s.l.) in PN-01 (Fig. 5a), An. barbirostris (s.l.) and An. dirus (s.l.) in PN-02 (Fig. 5b), An. dirus (s.l.) and An. donaldi in PN-03 (Fig. 5c) and An. barbirostris (s.l.) and An. dirus (s.l.) in PN-04 (Fig. 5d).

Plasmodium infection status
None of the 879 Anopheles mosquitoes tested were Plasmodium positive.

Anopheles mosquito parity status
A total of 3509 An. minimus (s.l.) and 322 An. baimaii were dissected to determine parity. Across the three districts the proportion of parous An. minimus (s.l.) and An. dirus (s.l.) ranged from 51.71 to 71.64% and from 55.32 to 74.71%, respectively (Table 3).

Power calculation
The purpose of this field study was to assess the potential of Surat Thani to serve as the location for future entomological evaluation of ivermectin MDAs. As Plasmodium infection rates in Anopheles in the GMS are extremely low, either survival or parity become the ideal metrics to measure. However, capturing indoor-resting, blood-fed Anopheles in the GMS is very difficult, rendering assessment of mosquito survival as done previously in Africa [12,13] impossible. The power calculation to assess the impact of ivermectin MDA on mosquito parity was based on An. minimus (s.l.) as it is a primary malaria vector, the most abundant Anopheles species captured in Surat Thani (87.16%; 5035/5777) and the most ivermectin-sensitive malaria vector in the GMS evaluated to date [16,17]. Parity rates (P 0 = 67.68%; 2375/3509) from ten clusters wherein > 100 An. minimus (s.l.) were collected and dissected were used to calculate an intra-cluster correlation (ICC) value to be used for sample size calculations. Individual mosquito parity results were assessed with a random effects logistic regression model to obtain an ICC of 0. 0378. To demonstrate an effect, based on a twosided α = 0.05, a parity rate at baseline and in control villages of P 0 = 67.68%, a treatment effect size = 34% reduction, meaning treatment villages with a parity rate of P 1 = 44.67% after MDA, and an ICC = 0.0378 requires at least four clusters per treatment arm and a minimum of 300 mosquitoes dissected per cluster to provide 82% power. A conservative treatment effect size of 34% was predicted based on a previously developed model for ivermectin MDA [41] parameterized with An. minimus (s.s.) ivermectin susceptibility data [17]. This treatment effect size is reasonable and supported by the following: (i) ivermectin MDAs in West An. dirus s.l.

Discussion
This study represents the most intensive and detailed analysis of Anopheles bionomics from Surat Thani province, Thailand. These surveys were conducted to determine which clusters in Surat Thani were most appropriate for evaluation of the entomological impacts of ivermectin MDA on Anopheles. Most (95.26%) of the Anopheles species collected from Phanom, Khiri Rat Nikhom and Vibhavadi districts were primary or secondary malaria vectors. In addition to detailed hourly collection data, this study presents data on the parity status throughout the night for the primary malaria vectors, An. dirus (s.l.) and An. minimus (s.l.). Finally, a power calculation was performed to determine the number of clusters that should be utilized to assess the impact of ivermectin MDA on An. minimus parity.
Surat Thani is unique in that it is a malarious province in Thailand that is not located along an international border. It has some of the highest rainforest coverage in the GMS and the highest rubber plantation coverage in Thailand, making it an ideal habitat for malaria vectors in the GMS. Indeed, of the 15 species captured in these surveys, 11 7.15% (407/5777) of the total Anopheles collected., respectively. Molecular identification confirmed that of all An. minimus (s.l.) collected, 99.80% were An. minimus (s.s.) (n = 484) and 0.2% were An. aconitus (n = 1), and of all An. dirus (s.l.) collected, 100% were An. baimaii (n = 348). Since the secondary and suspected vectors were captured in such low numbers, it is likely that the primary malaria vectors, An. minimus (s.s.) and An. baimaii, are largely responsible for Plasmodium transmission in the study area. However, as no specimens were Plasmodium positive (0/879), it cannot be confirmed that these vectors are solely responsible for transmission. A previous study in Ubon Ratchathani province, in northeastern Thailand, reported that suspected malaria vectors An. barbirostris (s.l.), An. philippinensis and An. hyrcanus (s.l.) are highly zoophagic, feeding mostly on cattle [43]. Since we only used the HLC method, this feeding preference may explain the low abundance of these species in the study area, making it difficult to rule them out as possible contributors to Plasmodium transmission in Surat Thani. One limitation of this study is the lack of assessment on the proximity of the mosquito collection sites to potential influential factors, such as the forest, larval habitat and livestock populations.
Primary malaria vectors in the GMS tend to be collected more frequently outdoors than indoors [9][10][11].
Since the aim of this study was to collect as many human host-seeking Anopheles as possible, all HLCs were performed outdoors. Results from a cross-sectional survey in Surat Thani in the same study area where HLCs were performed indicated that staying outdoors is a primary risk factor for asymptomatic Plasmodium carriage, suggesting that most of the transmission occurs outside the home [26]. This finding reinforces the potential usefulness of ivermectin MDA in the GMS as it can target the Anopheles malaria vector regardless of location or time.
Vibhavadi had the fewest number of Anopheles species collected, with only six species, and Anopheles species diversity was similar to that in Khiri Rat Nikhom, which were both lower than species diversity in Phanom (Table 1). This reduced Anopheles diversity in Vibhavadi could be due in part to seasonality and the limited sampling duration from July to October, while mosquitoes were collected from February to October in Khiri Rat Nikhom and Phanom. Another interesting point about mosquito collections in Vibhavadi was that An. minimus (s.l.) were captured later in the night compared to Phanom (P = 0.0275) and Khiri Rat Nikhom (P = 0.0002). This could be due to the production of durian in Vibhavadi as thermal fogging with malathion occurs at night-time for the control of several durian crop pests, in particular Scirtothrips dorsalis (Order Thysanoptera; Family Thripidae) and Allocaridara malayensis (Order Homoptera; Family Psyllidae), and this fogging was observed at some of the mosquito collection sites during HLCs in Vibhavadi clusters. Vibhavadi has the fourth highest coverage (1.33%) of durian plantation of the 19 districts in Surat Thani, with Khiri Rat Nikhom having the eighth highest coverage (0.33%) and Phanom the twelfth highest coverage (0.18%) [44].
The abundance of An. dirus (s.l.) varied among districts, with the greatest numbers captured in Phanom, likely due to geographical (Figs. 1, 2) and biological characteristics as the Phanom study site is surrounded by national parks comprised of primary old growth rain Table 4 Trend in parity by hour and cluster for Anopheles minimus (s.l.) (n = 3400) n number of mosquitoes dissected, OR odds ratio a All clusters showed an increasing trend in the proportion of parous An. minimus (s.l.) throughout the night (Fig. 10