High transmission efficiency of the simian malaria vectors and population expansion of their parasites Plasmodium cynomolgi and Plasmodium inui

Background The elimination of malaria in Southeast Asia has become more challenging as a result of rising knowlesi malaria cases. In addition, naturally occurring human infections with other zoonotic simian malaria caused by Plasmodium cynomolgi and Plasmodium inui adds another level of complexity in malaria elimination in this region. Unfortunately, data on vectors which are responsible for transmitting this zoonotic disease is very limited. Methodology/Principal findings We conducted longitudinal studies to investigate the entomological parameters of the simian malaria vectors and to examine the genetic diversity and evolutionary pattern of their simian Plasmodium. All the captured Anopheles mosquitoes were dissected to examine for the presence of oocysts, sporozoites and to determine the parous rate. Our study revealed that the Anopheles Leucosphyrus Group mosquitoes are highly potential competent vectors, as evidenced by their high rate of parity, survival and sporozoite infections in these mosquitoes. Thus, these mosquitoes represent a risk of human infection with zoonotic simian malaria in this region. Haplotype analysis on P. cynomolgi and P. inui, found in high prevalence in the Anopheles mosquitoes from this study, had shown close relationship between simian Plasmodium from the Anopheles mosquitoes with its vertebrate hosts. This directly signifies the ongoing transmission between the vector, macaques, and humans. Furthermore, population genetic analysis showed significant negative values which suggest that both Plasmodium species are undergoing population expansion. Conclusions/Significance With constant microevolutionary processes, there are potential for both P. inui and P. cynomolgi to emerge and spread as a major public health problem, following the similar trend of P. knowlesi. Therefore, concerted vector studies in other parts of Southeast Asia are warranted to better comprehend the transmission dynamics of this zoonotic simian malaria which eventually would aid in the implementation of effective control measures in a rapidly changing environment.


Introduction
Malaria is one of the most common vector-borne diseases in the world and it is endemic in many countries in tropical and subtropical regions [1]. Southeast Asia is in the pipeline for malaria elimination by 2030 [2]. Increasing cases of zoonotic simian malaria caused by P. knowlesi pose a new challenge to malaria elimination [3,4]. In addition, cases of natural human infection of other zoonotic simian malaria caused by P. cynomolgi [5][6][7][8][9][10][11] and P. inui [9,[12][13][14] have added another dimension of complexity to malaria elimination in Southeast Asia. Increasing cases of knowlesi malaria also hinder some of the Southeast Asian countries such as Malaysia from obtaining the malaria-free status certification from the World Health Organization (WHO), though the country had successfully eliminated indigenous human malaria cases since 2018 [15].
One of the contributing factors for the increasing cases of knowlesi malaria and possible emergence of other zoonotic simian malaria is due to extensive deforestation for agricultural activities and expansion of human settlements near forest fringes [16,17]. This potentially brought the spill over of the macaque population to human settlement which eventually caused the mosquito vectors to follow their macaque hosts and adapt to the new environment in semi-urban areas [18]. Unfortunately, the conventional measures to control human-malaria using long-lasting insecticide-impregnated bed nets (LLINs) and indoor residual spraying (IRS) are ineffective in the control of P. knowlesi transmission [19]. This is especially true since the main vectors for knowlesi malaria are Anopheles mosquitoes from the Leucosphyrus Group, which are known to be forest-dwelling mosquitoes [18]. Due to the exophagic and exophilic behaviour of the vectors, outdoor transmission and infective bites just after dusk and in the early morning remain a significant challenge for prevention and control [20][21][22]. Thus, understanding vector biology is paramount in strategizing effective control measures to combat the rising threat of zoonotic simian malaria. Explicitly, vector distribution, vector competency, adult behaviour and abundance play an important role in determining the transmission potential of the vectors [20].
Despite increasing cases of knowlesi malaria and possible emergence of other zoonotic simian malaria in Southeast Asia, there are limited vector studies to understand the transmission dynamic of this disease. Most of the vector studies have been focused on Malaysian Borneo [18], followed by Vietnam [23,24]. In Peninsular Malaysia, vector studies were conducted in selected locations a decade ago [25,26], and changes in the landscape due to extensive development for the past few decades may have altered the ecology and composition of the Anopheles mosquitoes. Thus, there is a dire need for an updated entomological study in this region given the drastic landscape changes due to substantial urbanization. Indeed, the impact of forest disturbance on changes in vector ecology is undeniable and it was well-studied in Malaysian Borneo [16].
In addition, comprehensive entomological studies are also warranted since most of the current vector studies relied heavily on PCR to detect the presence of Plasmodium parasites without mosquito dissection. As a result, crucial information such as mosquito parous rate is frequently omitted, which is the key index to estimate mosquito longevity and vector competency [27,28]. Sporozoite rate and entomological inoculation rate were also not deciphered [28]. Thus, there is a major knowledge gap especially on the entomological characteristics of the vectors in Southeast Asia. Furthermore, there is very limited information available regarding the prevalence and genetic diversity of the simian Plasmodium isolated from Anopheles mosquitoes.
Therefore, with increasing cases of zoonotic simian malaria, it is imperative to conduct a comprehensive entomological study on the local vectors and molecular studies on the simian malaria parasites. Thus, the present longitudinal study was conducted in Malaysia to determine the bionomics and transmission efficiency of the major vectors involved in the transmission of the zoonotic simian malaria parasites, and to evaluate the genetic diversity and evolutionary pattern of the simian Plasmodium in them.

Ethics statement
This study was approved by Medical Research and Ethics Committee, Ministry of Health Malaysia (NMRR-19-962-47606). Prior to the study, all volunteers who were involved in the mosquito collections had signed informed consent forms. The participants were asked to immediately report in case they felt ill or developed any symptoms. The project provided free blood examination and treatment of malaria for those who felt ill or wished to check themselves.

Study area
The study was conducted in seven different states in Peninsular Malaysia between June 2019 and March 2022 (Fig 1) 2"E) in Johor were chosen as longitudinal sampling locations. These locations were fixed for entomological investigations after preliminary study showed presence of Anopheles mosquitoes from the Leucosphyrus Group which were also positive for the presence of simian malaria parasites. Besides, other random sampling locations were selected based on past human cases of knowlesi malaria and through discussion with the district health officers (S1 Table). In total, there were 81 sampling locations, inclusive of the five longitudinal sampling locations throughout the study.

Mosquito collection
Adult female Anopheles mosquitoes were collected between June 2019 and November 2021 across all the sampling locations using human landing catch (HLC) and Mosquito Magnet as described in previous study [29]. For the longitudinal sampling locations, mosquitoes were collected for two consecutive nights per month from 1900 hours until 2300 hours for four months on a rotational basis between the longitudinal study locations. A total of 8 nights were spent at each longitudinal study location except for Bukit Tinggi, Johor where 10 nights were spent in mosquito collections by a team of two to three trained personnel. All collected https://doi.org/10.1371/journal.pntd.0011438.g001 Anopheles mosquitoes were carefully transported to the field laboratory for identification and were dissected to extract the midgut, salivary glands and ovaries to examine for oocysts, sporozoites and parity status respectively.

Mosquito identification
All the collected Anopheles mosquitoes were morphologically identified to the species level by using the taxonomic keys of Reid [30] and Sallum [31]. For Anopheles mosquitoes from the Leucosphyrus Group and mosquitoes which were difficult to be identified morphologically, they were confirmed molecularly through DNA sequencing by amplifying the ITS2 gene using primers ITS2A and ITS2B [32] with protocol as described in the previous study [33].

Statistical analysis for entomological indicators
All statistical analysis was carried out using R programming language (version 4.0.0). Parameters such as the abundance of mosquitoes, man-biting rates, parous rate and the proportion of mosquitoes which were infected with sporozoites and oocysts were analyzed using Generalized linear mixed model (GLMM). GLMMs were constructed in R using the glmmTMB package to compare the entomological parameters of An. introlatus from different sampling locations. In all the analyses, the localities were fitted as a fixed effect while the months of sampling were fitted as random effect. The entomological parameters for An. latens (n = 23) and An. cracens (n = 37) were not examined due to overall low density from the sampling locations for robust description.
Poisson distribution was assumed in the analysis of abundance and the man-biting rate of the mosquitoes. On the contrary, binomial distribution was assumed in the analysis of sporozoite rate, oocyst rate and the proportion of parous mosquitoes. Models testing associations between dependent variables (Vector abundance, man-biting rates, parous rates and infection rates) with fixed effect (sampling location) and random effects (months) were compared using higher log-likelihood and lower Akaike information criterion (AIC) values, along with the results of analysis of variance (ANOVA) of nested models. A Tukey's post-hoc test using multcomp package was used to assess the statistical differences of the entomological parameters between different localities.

Isolates used for 18S SSU rRNA gene characterisation
All the mosquitoes detected positive for the five simian Plasmodium were subjected to nested PCR targeting longer fragments of the 18S SSU rRNA gene for molecular characterisation. Since SSU rRNA gene is a widely used gene marker in eukaryotic phylogeny, there is an abundance of Plasmodium species sequence data available [34]. Thus, molecular characterisation of the SSU rRNA gene of the simian Plasmodium isolated from mosquitoes will allow comparison of sequence data from both human and macaques to better understand the transmission dynamic of the zoonotic simian malaria in this region. Besides mosquitoes, three macaques infected with simian Plasmodium (P. cynomolgi and P. inui) from each state were also included in this study. However, in Perak, only one macaque was detected to be positive for P. cynomolgi. The macaque samples were randomly selected from those captured closest to the mosquito longitudinal sampling locations. All the macaques were collected by the Department of Wildlife and National Parks of Peninsular Malaysia (PERHILITAN) and screened for the presence of simian Plasmodium in previous study [35]. In that study, both, P. inui and P. cynomolgi were discovered in high prevalence among the macaques in Peninsular Malaysia. The positive samples were then subjected to nested PCR targeting longer fragments in this current study for molecular characterisation together with the DNA sequences of simian Plasmodium amplified from mosquitoes.

Detection of Plasmodium parasites in infected mosquitoes and nested PCR assay for 18S SSU rRNA gene characterisation
The mosquitoes were examined for the presence of sporozoites in the salivary glands and for oocysts in the midgut. Genomic DNA was extracted from the parasite-positive guts and glands using the DNeasy tissue kit (Qiagen, Germany) according to the manufacturer's protocol. Nested PCR assay was performed targeting the Plasmodium small subunit ribosomal RNA (18S rRNA) gene to identify human malaria parasites (Plasmodium falciparum, P. malariae, P. ovale curtisi, P. ovale wallikeri and P. vivax) and simian Plasmodium (P. coatneyi, P. cynomolgi, P. fieldi, P. inui and P. knowlesi) using genus-specific primers rPLU 1 and rPLU 5 for the nest 1 amplification [36], followed by species-specific primers in the nest 2 amplification [37][38][39] (S2 Table). This protocol has been described in a previous study [29]. The amplification products were analyzed using 1.5% agarose gel electrophoresis.
For gene characterisation targeting the large fragment of the 18S SSU rRNA gene, the same nest 1 product was used. PCR amplification reaction for nest 2 assay was performed using a universal forward primer UMSF [40] combined with species specific reverse primers [38][39][40]. PCR amplification reaction for nest 1 assay was performed in a final volume of 50 μL containing 5 μL of DNA template, 1× GoTaq Long PCR Master Mix (Promega, USA), 0.25 μM each of forward (rPLU1) and reverse (rPLU5) primers. Cycling parameter for nest 1 consisted of initial denaturation at 94˚C for 4 min, followed by 40 cycles of 94˚C for 30s, 55˚C for 1 min, 72˚C for 1.5 min and a final extension at 72˚C for 10 min. For each 25 μL of nest 2 PCR amplification, 3 μL of nest 1 PCR amplification product was used as DNA template. PCR amplification reaction for nest 2 assay consisted of 1× Green GoTaq reaction buffer (Promega, USA), 2.0 mM MgCl 2 (Promega, USA), 0.2 mM of dNTPs mixture (Promega, USA), 0.4 μM each of forward and reverse primers and 1.5 U of GoTaq DNA polymerase (Promega, USA). Cycling parameter for nest 2 consisted of initial denaturation at 94˚C for 5 min, followed by 35 cycles of 94˚C for 1 min, species specific annealing temperature for 1 min, 72˚C for 1 min and a final extension at 72˚C for 5 min (S3 Table).

Cloning and sequencing of Plasmodium 18S SSU rRNA gene fragments
The large fragment of 18S SSU rRNA gene of the simian Plasmodium species isolated from the Anopheles mosquitoes and macaques were cloned and sequenced. The amplified PCR product was excised from the agarose gel using a sterile scalpel after gel electrophoresis. The gel slice was purified using NucleoSpin Gel and PCR Clean-up (Macherey-Nagel, Germany) according to the manufacturer's protocol. Cloning the SSU rRNA gene was done using pEASY-T5 Zero cloning kit (TransGen Biotech, China) according to the manufacturer's protocol. For transformation, Trans1-T1 Phage Resistant Chemically Competent Cell (TransGen Biotech, China) was added to the ligated product. After incubation period, colony PCR was conducted to identify the positive clones with gene insert. Colony PCR was conducted using M13F (F-20) forward primer and M13R (R-26) reverse primer. After the propagation step, the plasmid DNA was extracted using QIAprep Spin Miniprep kit (Qiagen, Germany) according to the manufacturer's protocol. The extracted plasmids were then sent to 1 st Base Laboratories Sdn. Bhd, Malaysia for DNA sequencing.

Sequence editing, alignment and phylogenetic analysis
The nucleotides sequences of SSU rRNA gene of the simian Plasmodium parasites obtained from this study were aligned and trimmed using BioEdit version 7.2 software (https://bioedit. software.informer.com/7.2/) while similarity searches were conducted using Basic Local Alignment Search Tool (BLAST) (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Nucleotide sequences of the SSU rRNA gene obtained from this study were phylogenetically compared with those obtained from GenBank for the five simian Plasmodium species using MEGA version 10.1 software (https://www.megasoftware.net/). Phylogenetic tree was constructed using neighbour-joining method and the evolutionary distances were computed using maximum composite likelihood model with a bootstrap value of 1000 replicates to test the robustness of the tree [41]. In the phylogenetic tree, Plasmodium berghei (AJ243513.1) was used as an outgroup.
Nucleotide sequences obtained from cloned samples were aligned with ClustalW using MEGA version 10.1 software (https://www.megasoftware.net/). Sequence analysis and nucleotide comparison were carried out against reference sequence for P. inui San Antonio I strain (GenBank accession number: XR606809) and P. cynomolgi Mulligan strain (GenBank accession number: AB287290). For P. inui and P. cynomolgi, 171 sequences of PinA-type 18S rRNA (933 bp) and 73 sequences of PcyA-type 18S rRNA (927 bp) were respectively used for analysis. The accession number of imported DNA sequences were listed in S4

Haplotype network analysis
The DNA polymorphism of the A-type 18S SSU rRNA genes for P. cynomolgi and P. inui were estimated by calculating the number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), number of polymorphic sites (k) and average number of pairwise nucleotide differences using DnaSP version 6.12.03 which allows comprehensive analysis of DNA sequence variation. Haplotype networks were also constructed for both Plasmodium species based on the A-type 18S SSU rRNA partial sequences generated from the present study as well as published sequences from GenBank databases by using the median-joining method in NET-WORK version 10.2.0.0 software (Fluxus Technology Ltd, UK). Where available, DNA sequences of reference strains of P. cynomolgi and P. inui were included in the construction of the haplotype networks. Details of the published sequences extracted from GenBank databases used in the analysis were listed in S4 and S5 Tables for P. cynomolgi and P. inui respectively.

Population genetic analysis
The 18S SSU rRNA gene was analysed to elucidate the population expansion of simian Plasmodium isolated from Anopheles mosquitoes. Populations pairwise F ST values for genetic distance between the subpopulations were tested for significance using DnaSP version 6.12.03. The value is ranged from 0 to 1 and interpreted as F ST = 0 (no differentiation), F ST < 0.05 (little differentiation), F ST between 0.05 to 0.15 (moderate differentiation), F ST between 0.15 to 0.25 (great differentiation) while F ST > 0.25 (very great differentiation) [42]. The gene flow was categorized as Nm < 0.25 (low gene flow), Nm between 0.25 to 0.99 (intermediate gene flow) while Nm > 1 (high gene flow) [43]. On the other hand, neutrality test was conducted using Tajima D test [44], Fu and Li D* [45] and Fu and Li F* [46] statistics by using DnaSP version 6.12.03. Under neutrality, Tajima's D is expected to be 0. Significantly positive Tajima's D values indicate recent population contraction with selection maintaining variation, whereas negative values suggest population expansion with selection removing variation [44]. Demographic expansion was further investigated with mismatch analysis test using raggedness index (Rag) [47].

Abundance of Anopheles species
A total of 1652 Anopheles mosquitoes belonging to 15 species were caught and examined for the presence of both human and simian Plasmodium parasites (Table 1). Anopheles maculatus was the predominant species obtained comprising 41.9% of the total catch, followed by An. introlatus (24.3%), An. letifer (12.7%) and An. sinensis (11.1%). Other Anopheles species only recorded less than 10% of the total catch. Among the 15 Anopheles species obtained throughout the study, An. cracens, An. introlatus and An. latens were the only mosquito species belonging to the Leucosphyrus Group.

Species composition of Anopheles in all longitudinal study locations
A total of 607 Anopheles mosquitoes belonging to 11 different species were collected from five longitudinal study locations across Peninsular Malaysia ( Table 2). From those longitudinal sampling locations, An. introlatus was the predominant species collected which comprised 36.4% of the total catch, followed by An. letifer and An. maculatus, each representing 25.9% of the total catch. Other Anopheles species were found in lower abundance where each species recorded less than 10% of the total catch.

Seasonal abundance and biting rate of Anopheles from longitudinal sites
The number and species of Anopheles mosquitoes vary across different sampling months in each longitudinal sampling location (Fig 2). By using the GLMM models, the results demonstrated that the Poisson distribution was generally a better model than the binomials and negative binomials in analyzing the abundance and biting rate of An. introlatus. Tukey post hoc test revealed that the abundance and biting rate of An. introlatus were significantly higher in Bukit Tinggi (forest) compared to other longitudinal sampling locations (P<0.05). By controlling the effect of monthly variation, the mean predicted abundance of An. introlatus in Bukit Tinggi (forest) was 19 to 28 times higher than other sampling locations, with 12 to 20 times higher mean predicted biting rate. An. barbirostris gp. Different peak biting times were observed for the different Anopheles species collected from this study (Fig 3). The peak biting time of An. introlatus was between 2000 to 2100 hours followed by a sharp decline in the percentage of mosquito biting thereafter. Interestingly, the same trend was observed for all An. introlatus from different longitudinal sampling locations. For An. latens and An. cracens, the peak biting time was observed slightly earlier than An. introlatus which was between 1900 to 2000 hours. This was observed from Kg. Lalang (village) and Gunung Panti (forest) where An. latens were found and in Kem Sri Gading (forest) for An. cracens. Generally, all Anopheles mosquitoes from the Leucosphyrus Group from this study were early biters with high biting rate between 1900 to 2100 hours with a steady decline thereafter.

Transmission efficiency of the Anopheles Leucosphyrus Group
The parous rates for all three species of Anopheles mosquitoes from the Leucosphyrus Group (An. cracens, An. introlatus and An. latens) were predominantly more than 60%. Nonetheless, seasonal variation in the parous rate was not observed for all three species in longitudinal sampling locations (Fig 4). The mean parous rate for An. introlatus ranged from 68.32% to 75.00% which was slightly higher than An. latens (66.67% to 70.00%) and An. cracens (66.67%)   (Table 3). Tukey's test showed no significant difference in the parous rate of An. introlatus collected from different sampling locations. Statistical test was not employed for the other two species of Anopheles since they were only found in some of the longitudinal sampling locations, thus precludes any meaningful comparison. Estimation on the probability of daily survival rate, life expectancy and vectorial capacity based on the parous rate were computed for all Anopheles mosquitoes from the Leucosphyrus Group at each sampling location (Table 3). Generally, the life expectancy of An. introlatus was higher compared to An. latens and An. cracens from this study. It would be expected that between 28% to 38% of An. introlatus would survive the 10 days necessary for the formation of P. knowlesi or other simian Plasmodium sporozoites (7.5 days for P. cynomolgi and 11 days for P. inui) [48]. Those surviving the 10 days would have a further infective life expectancy of 2.2 to 4.0 days ( Table 3). The vectorial capacity was the highest for An. introlatus from Bukit Tinggi (forest) (4.60) compared to An. introlatus from other longitudinal sampling locations, mainly influenced by the high man biting rate. The low number of An. introlatus collected from other longitudinal study locations besides Bukit Tinggi also partially influenced the vectorial capacity of the mosquito.

Infection rates and entomological inoculation rates by months and sites
The oocyst rate and sporozoite rate of the three known simian Plasmodium vector species (An. cracens, An. introlatus and An. latens) were assessed. There was no consistent seasonal pattern observed for the infection rates for all the three Anopheles species across different study locations (Fig 5). For An. introlatus, the oocyst rate was extremely high in the month of December 2020 (100%) at Sungai Dara (forest). This was mainly due to low sample size because only a single An. introlatus was captured and dissected and was found to be positive for the presence of oocysts. A similar scenario was also observed for An. latens in the month of August 2019 in Kg. Lalang (village) where the oocyst rate was recorded 100%. GLMM analysis with Tukey post hoc tests on An. introlatus infection rates indicated that there was no significant difference between the infection rate of An. introlatus among the five different sampling locations. Limited An. latens (n = 23) and An. cracens (n = 37) were collected for robust statistical analysis. These two Anopheles species were only found in some of the longitudinal sampling locations, thus making it impossible for comparative analysis between different study locations. Overall, the collected Anopheles mosquitoes from the Leucosphyrus Group had an entomological inoculation rate (EIR) between 0.04 to 0.27 (Table 3 and Fig 5). The high EIR was observed in Bukit Tinggi (forest) (0.27) mainly due to high man biting rate. This was followed by An. cracens at Kem Sri Gading (forest) (0.12) while the lowest EIR was recorded by An. introlatus from Gunung Panti (forest) (0.04).

Identification of Plasmodium species in mosquito samples
All the collected Anopheles species (n = 1652) were examined for the presence of Plasmodium parasites through dissection of the mosquitoes to detect the presence of sporozoites or/and oocysts. For those mosquitoes which could not be dissected, the head and thorax were separated from the abdomen, and they were individually screened using PCR to detect the presence of Plasmodium parasites. All the positive mosquitoes were mainly from the Leucosphyrus Group except for An. letifer which was from the Umbrosus Group.
The predominant Plasmodium species recovered in this study were P. inui (n = 33), followed by P. fieldi (n = 10), P. cynomolgi (n = 9) and P. coatneyi (n = 2) ( Table 4). For the successfully identified Plasmodium species, basic local alignment search tool (BLAST) had 97-100% similarities between all the sequences acquired from this study against the sequences available in GenBank. Most of the positive mosquitoes were of mono-infection (80%), followed by double (18%) and triple infections (2%). The single infections were mainly P. inui (n = 23). However, the parasite from some mosquito samples, An. letifer (n = 4), An. introlatus (n = 2) and An. latens (n = 1) could only be identified to Plasmodium and were negative for human and the five simian malaria parasites species.

Phylogenetic analysis
Phylogenetic tree constructed based on 18S SSU rRNA gene of the simian Plasmodium isolated from Anopheles mosquitoes revealed five monophyletic clades with bootstrap values ranging from 59% to 100% (S1 Fig). All the 103 simian Plasmodium clones from 54 Plasmodium isolates identified from this study were assigned in the phylogenetic trees with their respective species group. The five major clades formed in the phylogenetic tree corresponds to the five simian Plasmodium species commonly identified in Malaysia, namely P. coatneyi, P. cynomolgi, P. fieldi, P. inui and P. knowlesi. Besides, some of the previously published sequences of simian Plasmodium from Malaysian Borneo were also included in the phylogenetic tree (26 sequences). These samples were positioned in the same clade together with samples originating from Peninsular Malaysia based on the species group.

Haplotype network analysis
The haplotype network of P. cynomolgi (PcyA-type 18S rRNA) gene showed geographical clustering between Peninsular Malaysia and Malaysian Borneo which formed two dominant haplotypes, H_8 and H_4 respectively (Fig 6A). The haplotype H_8 was shared between the P.

Population genetic structure
Overall, the neutrality tests for both P. cynomolgi and P. inui populations showed significant negative values which indicate that the populations are undergoing expansion (S6 Table). In contrast, neutrality tests according to different hosts for P. cynomolgi isolates from Peninsular Malaysia had shown positive values suggesting recent population contraction due to balancing selection. However, the values were not statistically significant except for Tajima D test for P. cynomolgi isolated from mosquitoes in Peninsular Malaysia (P < 0.05). For P. inui, neutrality tests on isolates of Anopheles mosquitoes from Peninsular Malaysia showed significant negative values suggesting recent population expansion (P < 0.01). The pairwise distribution analysis for both P. cynomolgi and P. inui isolates from Peninsular Malaysia revealed unimodal distribution suggesting population expansion (Fig 7). This was in contrast with the neutrality  Table). Similarly, P. inui populations isolated from mosquitoes in Peninsular Malaysia showed a very great genetic differentiation compared to macaques (F ST = 0.683), humans (F ST = 0.683) and mosquitoes (F ST = 0.601) from Malaysian Borneo. Besides, a high F ST values were also obtained when comparing Malaysian populations with populations from other countries such as Thailand (F ST > 0.25) and Taiwan (F ST > 0.40), high likely due to the geographical distance which limits the gene flow between the subpopulation of P. inui (S8 Table).

Discussion
With increasing cases of zoonotic simian malaria in Southeast Asia, it is of paramount importance to identify potential vector species and to understand their bionomics to elucidate human exposure risk to this zoonotic simian malaria. In this study, An. introlatus was the predominant Anopheles species (36.4%) collected from the longitudinal surveillance with mean biting rates ranging from 0.33 to 6.30 bites/man/night. Anopheles mosquitoes belonging to the Leucosphyrus Group were the predominant group in each longitudinal sampling location except for Gunung Panti (forest) where the Umbrosus Group was the dominant group. The predominant Anopheles species varies across different geographical locations not only in Peninsular Malaysia but also in Southeast Asia. Variation in the species composition and abundance can perhaps be attributed to the environmental heterogeneities [49]. This could also explain why the cases of simian malaria vary widely throughout the region. A preliminary study conducted in Bukit Tinggi army camp in the forest revealed that none of the Anopheles mosquitoes were caught indoors while significantly high number of Anopheles mosquitoes were successfully collected outside the camp. This is in parallel with previous studies which showed the exophagic behaviours of Anopheles mosquitoes from the Leucosphyrus Group [26,50]. Thus, indoor collection was not carried out throughout the duration of this study. A substantial proportion of Anopheles mosquitoes from this study were found biting outdoors in the early evening from 1900-2100 hours, a time when most inhabitants in villages are actively engaged in outdoor activities and unprotected by insecticides treated bed nets [51]. There is also a higher chance of human being exposed to infectious bites when they are returning from daily activities in the forest [52].
Generally, the parous rate for Anopheles mosquitoes from this study was high, ranging from 66.67% to 75.00%. The high percentage of parous rate might possibly be due to the suitable climate conditions such as high relative humidity in the area which has a significant impact on the mosquito longevity [27], especially when all the sampling locations were either in forested area or forest fringes with dipterocarp trees forming canopies which provides shade and suitable humid environment for the mosquitoes to thrive. Calculating the parous rate is an important parameter that reflects the age or life expectancy of the mosquito population which has its epidemiological importance in malaria disease transmission [53]. It is only the parous mosquitoes that are likely to transmit malaria because it provides sufficient time for the parasites to develop until the infective sporozoite stage [54].
Although other simian malaria is being reported in small numbers in humans, it should not be dismissed because natural transmission of these zoonotic simian Plasmodium is already starting to take root. This includes P. inui infection in human in Malaysia [13] and P. fieldi infection in Thailand [12]. While in Kalimantan Indonesia Borneo, human cases of simian malaria caused by P. coatneyi, P. fieldi and P. inui have been postulated [14], which indeed highlight the potential of these simian Plasmodium to infect humans in natural settings.
In general, this study revealed that P. inui was the predominant Plasmodium species followed by P. fieldi and P. cynomolgi respectively. This observation is in parallel with the prevalence of simian Plasmodium recovered from macaques in Malaysia [4]. The over-lapping distribution of the vectors and the macaques in forested areas might explain the high prevalence of these simian Plasmodium especially P. cynomolgi and P. inui in both the macaques and in the Anopheles mosquitoes. Unfortunately, P. knowlesi was not recovered in any of the Anopheles mosquitoes in this current study. Surprisingly, this is also in parallel with few studies in Sabah [55,56] and Sarawak [57] where extremely low number of P. knowlesi were detected in mosquitoes.
Nevertheless, failure to detect P. knowlesi in mosquitoes in this study should not be construed as evidence of no transmission. The inability to detect P. knowlesi in infected Anopheles mosquitoes from this study might be due to its low prevalence in mosquitoes as shown in many previous studies [4] or due to selection of sampling locations. Although the sampling locations were selected based on the past knowlesi malaria cases, most of the mosquito samplings were conducted at forest fringes or forested areas accessible by vehicle. However, the infections might have occurred much deeper inside the forest since the villagers who are engaged in forest-related occupation have been reported staying overnight in the deep forest which are usually not accessible by vehicle. Indeed, studies have shown that forest related activities are one of the key risk factors for acquiring P. knowlesi infections [58]. In addition, a much larger sample size with broader sampling locations covering a range of forest sites may be required to accurately estimate the prevalence of P. knowlesi in Anopheles populations in Peninsular Malaysia. The ability of Anopheles vector to transmit other simian Plasmodium such as P. cynomolgi and P. inui besides P. knowlesi poses some threat through the emergence of these parasites in future. This is more likely since there were studies highlighting co-infections of these simian malaria parasites with human Plasmodium in mosquitoes [23,24,59] indicating possible simultaneous transmission. Since a substantial number of the vectors were positive for P. inui and P. cynomolgi, and asymptomatic cases are known to occur [6] it is crucial to screen people in such locations to ensure that they are free of the parasites.
On the other hand, the constructed haplotype network provided insight into the phylogeography clustering of the parasite populations [60]. The shared dominant haplotype of P. inui, H_18 between mosquito, macaques and human suggests evidence of parasite transmission between the hosts and human. This was further supported by the detection of two natural P. inui human infections in Pahang, Malaysia together with P. inui positive An. cracens mosquito from the same study site [13]. Similar scenario was also observed in the haplotype network of P. cynomolgi where human and macaques from Peninsular Malaysia shared the same haplotype H_8 while macaques and Anopheles from Borneo Malaysia shared another haplotype H_4. Clustering of haplotype of infected Anopheles mosquitoes along with the shared haplotypes between different hosts indicate the possible transmission between mosquitoes, humans and macaques in this region. This is especially true when the region had undergone extensive ecological changes due to deforestation giving way for oil palm plantations over the years [61]. Environmental changes especially associated with deforestation and land exploration bring human population in close contact with macaques and forest dwelling Anopheles mosquitoes over time, which inevitably increases the risk [62]. This is further supported by recent study on the blood meal of the vectors from the Leucosphyrus Group in this region which exhibit simio-anthropophagic behaviours [63].
Haplotype of P. inui from Taiwan and Thailand showed distinct geographical clustering away from haplotypes of P. inui isolated from Peninsular Malaysia and Borneo Malaysia. Indeed, P. inui isolates from Taiwan [64] and Thailand [65] were collected from a few macaques at sampling sites near each other in their respective countries which eventually leads to the formation of a tight clusters in the haplotype network. Contrarily, the DNA sequences used in the analysis from both Peninsular Malaysia and Borneo Malaysia were collated from the current study and a few different studies conducted many years apart, thus possibly forming a broader cluster [9,[66][67][68]. Besides, the separate clustering of haplotypes derived from macaques from Taiwan could also be due to different species of host where the parasite was isolated. The P. inui from Taiwan was isolated from Formosan macaques (Macaca cyclopis) [64] which was different compared to P. inui isolated from long-tailed macaques (Macaca fascicularis) from other countries. Indeed, a previous study on P. knowlesi using multilocus microsatellite genotyping had shown two divergent clusters associated with different species of macaque hosts (M. fascicularis and M. nemestrina) [69]. Although elucidative and in parallel with previous studies, the conclusion gained from the haplotype analyses remain speculative as they were derived from limited data especially for P. inui and P. cynomolgi isolated from humans and macaques. On the other hand, most of the simian Plasmodium isolated from vectors in this study originated from southern Peninsular Malaysia which had relatively high sample size of Anopheles mosquitoes. Thus, additional studies are required using larger sample size of these zoonotic simian Plasmodium species isolated from all the three different hosts from various geographical locations. It is also worth analyzing other polymorphic genetic markers especially from rapidly evolving genes such as Var gene [70] or microsatellite loci of the parasites [69] to observe a clear geographical clustering of the parasite populations.
The overall neutrality results for P. cynomolgi and P. inui revealed negative values, suggesting both the populations might be undergoing demographic expansion; a similar trend observed for both P. knowlesi [71] and other human Plasmodium species [72,73]. The unimodal mismatch distribution with low values of Raggedness index (P < 0.05) further supported the hypothesis of population expansion for both the simian Plasmodium species from this study.
With constant microevolutionary processes, these simian Plasmodium will soon be able to adapt and evade the immune responses in wild vector populations and perhaps alter the Plasmodium virulence in natural population. Indeed, the effects of genetic variation and changing environmental factors do play a vital role in altering malaria parasite infectivity and Anopheles susceptibility to these parasites which eventually affects the emerging infection rates in the future [74]. History has proven that non-human primates do share malaria parasites with humans and host-switching is an ongoing evolutionary process [75]. Thus, if "host-switching" events were to take place, P. knowlesi and other emerging zoonotic simian malaria caused by P. cynomolgi and P. inui might increase in their transmissions and eventually pose a greater threat to the public health in the future. Therefore, besides entomological studies on the vectors, more genetic study on these zoonotic simian Plasmodium are warranted to further bolster the understanding on the evolution of these parasites and its interactions between the hosts and vectors which affects its transmission capabilities. Since there is already evidence of simian malarias other than P. knowlesi being transmitted to humans, proactive measures are needed to prevent escalated risk of infection in the future especially when human malaria is eliminated.