Aedes aegypti is an efficient vector of dengue due to its highly adaptive nature to the urban environment. Although it is observed to have a short dispersal (active) capability, it has been shown to be capable of traveling long distances (passive) via human-mediated transportation. This duality may expand the distribution of the mosquito vector in urbanized areas. In this study, we examined the population genetic structure of Ae. aegypti in a highly urbanized and dengue-endemic region of the Philippines, Metropolitan Manila. Our findings indicated the dual dispersal nature of Ae. aegypti. The use of microsatellites as genetic markers also allowed us to describe the potential long-distance dispersal patterns, possibly through human-aided land transportation via the existing road networks of Metropolitan Manila.
Dengue disease is the most prevalent mosquito-borne viral infection in tropical and subtropical countries  with approximately 2.5 billion people at risk of contracting the disease worldwide . The Philippines considers dengue as a major public health concern and has been ranked fourth in Southeast Asia with the highest number of dengue cases reported annually . In particular, the National Capital Region (NCR), also known as Metropolitan Manila, is one of the main areas affected by the disease and an average of 16% of the total number of cases reported in the Philippines from 2009 to 2014 . Interestingly, the first recorded dengue epidemic in Southeast Asia occurred in this region in 1954 . Although a dengue vaccine is available , the World Health Organization  still recommends disease prevention and control towards the mosquito vector.
The transmission of this arboviral disease is through the bite of either the Aedes aegypti or Ae. albopictus mosquito vector. The former is identified as the primary vector of the disease, and as such, its biology has been extensively studied including its reproduction and development [7–9], survivability [10–11], dispersal [12–14], and ecology[15–16]. In many population genetic studies of Ae. aegypti , microsatellites have been used as a molecular marker of choice for differentiating macro- [17,19] and micro-geographic [20,21,22] scale mosquito populations due to its high polymorphism, co-dominance, and broad genome distribution . Molecular genotyping of the mosquito vector using these markers has shown additional insights towards the microevolution of the mosquito vector [17–19], revealing the population’s gene flow pattern which can be interpreted as the mosquito vector’s dispersal ability [20,24].
Urbanization has been a key factor in the widespread occurrence of dengue disease worldwide and has contributed to the dispersal of the mosquito vector . Ae. aegypti has been documented to thrive in highly urbanized areas, creating the ideal condition for a dengue epidemic [25,26]. As such, investigating the impact of urbanization on Ae. aegypti ’s population genetic structure has sparked research interest, particularly regarding its micro-evolutionary process [24, 26]. A growing urban area and human population density can affect the mosquito populations on a genetic level. An increase in urbanization can reduce genetic diversity and promote homogenization of the genetic structure of mosquito populations. However, when an area reaches complete and/or stable urbanization, the mosquito populations may reach peak abundance, allowing an increase in genetic diversity and thereby, establishing heterogenous genetic structures of mosquito populations ; a phenomenon that has been observed in urbanized cities in Brazil , Thailand , Indonesia , Cambodia and Myanmar .
Although dengue is endemic in the Philippines, there are a limited number of studies that investigate the population genetics of Ae. aegypti in the country. These studies revealed the seasonal genetic structuring and expansion of the mosquito vector  and the connectivity among its populations from different island seaports , but provide little insight into genetic structuring of Ae. aegypti in highly urbanized and dengue-endemic regions in the country. In this study, we selected Metropolitan Manila because it represents an ideal setting to investigate how urbanization can shape the present genetic structure of Ae. aegypti . Previous studies conducted in this area have only investigated the association between dengue epidemiology and climate [5,30,31] and the population structure using wing geometric morphometrics , offering little information about its genetic structure and relatedness. Therefore, the study collected adult Ae. aegypti samples from several sites within Metropolitan Manila in order to characterize the population genetic structure and relatedness among mosquito populations.
Metropolitan Manila, a highly urbanized area, is the National Capital Region (NCR) of the Philippines with an area of 636 km2 . It is located at the eastern shore of Manila Bay in Southwestern Luzon (14°50′ N Latitude, 121°E Longitude). It is composed of 16 cities and 1 municipality with a total population of 12,877,253 . This area is the most urbanized region in the Philippines and is the center of the national government, economy, education and culture, business, manufacturing and importation and exportation .
In this study, Metropolitan Manila was divided into 21 study areas (Fig 1 and S1 Table). Initially, study areas were delineated per city with at least one study area per city. However, some cities comprised either a large or small land sized area, thus designating more study areas. For example, the largest city, in the region Quezon City was divided into five study areas based on its district boundaries. In order to standardize the land size covered by each study area, we merged two small neighboring cities, San Juan and Mandaluyong. All maps were prepared in ArcGIS software version 10.2.2 (ESRI, NC, USA) for further analysis .
Households on each study area were selected based on voluntary informed consent for collecting adult Ae. aegypti mosquitoes inside their premises. The number of households per study area ranged from 2–14 with an average of six households in each area and a total of 151 households. The maximum distance between households within study areas ranged from 1.39–6.17 km. Since the households were widely dispersed across each study area, we calculated the geographical midpoint  to assign a single georeferenced location for each study area in the subsequent genetic analysis. The distance among study areas (midpoints) ranged from 2.85–39.66 km.
The majority of the previous population genetic studies of Ae. aegypti have either used only larval or reared larval-to-adult samples which could lead to potential bias in estimating population genetic parameters due to the sampling of full sibling mosquito larvae . This was evidenced in a population genetic structure study of amphibian larval samples that led to inaccurate estimates of differentiation among populations when compared to adult samples [38,39]. As such, this study targets adult Ae. aegypti samples rather than conventional egg, or larval samples to prevent collecting mosquito full siblings. Collection of Ae. aegypti mosquitoes was done by installing a commercially available mosquito UV Light Trap (MosquitoTrap, Jocanima Corporation, Las Pinas City, Philippines) in each household’s outdoor premises for 3–5 days. Collected adult mosquito individuals were sorted, then identified accordingly based on the pictorial keys of Rueda et al.  and preserved in 99% ethanol. From May 2014 until January 2015, a total of 526 adult Ae. aegypti were collected, ranging from 10 to 42 individuals per study area (S1 Table).
The total genomic DNA of individual mosquito sample was extracted using the QIAGEN Blood and Tissue DNEasy Kit (Qiagen, Hilden, Germany) following the modified protocol of Crane . The average DNA concentration per individual mosquito sample was measured at 41.68 ± 34.50 ng/ μL using NanoDrop2000 Spectrophotometer (Thermo Scientific Inc., Waltham, Massachusetts, USA). We selected 11 microsatellite markers out of the 24 loci as listed by Slotman et al.  and Chambers et al.  by conducting trial runs (PCR and gel electrophoresis) on a subset of Ae. aegypti samples (n = 8). The selection of the markers was based on the reliability of the amplified products in all samples from different trial runs. The selected markers used in this study were congruent with previous population studies conducted in Southeast Asia [19,21]. The loci were grouped into four sets for multiplex PCR (S2 Table). Each set consisted of fluorescent labeled forward primers with different annealing temperatures during the amplification process. Each set was composed of 1.2 μL of 10X buffer (TaKaRa, Shiga, Japan), 0.8 μL of 25 mM MgCl2 , 1.6 μL of 10 mM of each dNTP, 0.6 μL of 10 μM forward and reverse primers and 0.08 μL 5.0U/ μL of Taq DNA polymerase (TaKaRa), 1.5 μL of 10% dimethyl sulfoxide (DMSO) (Sigma-Aldrich, St. Louis, Missouri, USA) and 1 μL of template DNA with a final volume of 10 μL. Thermocycling conditions were as follows: initial denaturation at 94°C for 5 minutes, denaturation at 94°C for 30 seconds, annealing (see S2 Table for temperatures and durations for each primer set) extension at 72°C for 30 seconds (35 cycles) and a final incubation step of 72°C for 5 minutes. PCR amplicons were checked by electrophoresis in 3% agarose gels and visualized under UV light using the Chemidoc XRS Chemiluminescent Gel Documentation Cabinet (BIO-RAD). Prior to fragment size analyses, multiplex PCR products were diluted in 1/15 water and then pooled together. One microliter of each diluted pool was added with 0.5 μL of GS 500 Liz Internal Size Standard (Applied Biosystems, USA) and HD formamide for a total volume of 20 μL. Fragment analysis of the amplified products were done using ABI 3500 Genetic Analyzer (Life Technologies) while genotyping was conducted in GeneMapper (Applied Biosystems). Microsatellite data were checked for error, and the presence of null alleles was confirmed with MicroChecker . The dataset was deposited and made publicly available at vectorbase.org  with the population biology project ID: VBP0000554.
The exact Hardy-Weinberg equilibrium (HWE) test and estimations of the linkage disequilibrium (LD) among all pairs of loci were conducted using GENEPOP v4.2.1 [45,46]. Significance levels for multiple testing were corrected using the Bonferroni procedure. Measures of genetic diversity, including mean number of different alleles (MNa), mean number of effective alleles (MNe), allelic richness (AR), private allelic richness (PAR), observed heterozygosity (Ho), expected heterozygosity (He ), and inbreeding coefficients (Fis) were estimated using Genetic Analysis in Excel (GenAlEx) version 6.3 . To assess the magnitude of genetic differentiation among sites, uncorrected and corrected pairwise FST values were calculated using Arlequin v184.108.40.206  and FreeNA , respectively, with 10,000 permutations. The analysis of molecular variance (AMOVA) was performed using GenAlEx version 6.3 . Since the corrected pairwise FST values were calculated by FreeNA , we manually calculated the number of pairwise effective migrants using the formula of Slatkin and Barton  in GenAlEx version 6.3 . FreeNA  was also used to calculate the Cavalli Sforza and Edwards distance. A dendrogram was constructed based on the Cavalli Sforza and Edwards distance using the unweighted pair group with arithmetic mean (UPGMA) in the fastcluster package  and the optimal number of clusters were determined using the cindex index in the NbClust package  of the R program .
The number of genetic groups (K) was inferred using the Bayesian approach in the software Structure version 2.3.2 . The study used admixture models where the alpha value may vary, and independent allele frequencies set at lambda equals to one. Twenty independent runs were performed for each value of K (1–20) with a burn-in phase of 200,000 iterations followed by 600,000 replications. To determine the most likely number of clusters, we calculated the commonly used ΔK statistic as developed by Evanno et al (2005)  using the web application Structure Harvester version 0.6.93 . To avoid the effects of uneven sampling [57,58] in the Bayesian analysis using the software Structure version 2.3.2 , we standardized the number of individuals per study area to 10 and created five different datasets. These datasets were subjected to Bayesian analysis using the same parameters and analyzed using the Evanno et al (2005)  method. The ‘pophelper’ package  in R  was used to summarize and visualize the bar plots for the best K statistic identified for each dataset.
Pairwise geographic distances (km) among study areas and households of the mosquito populations were calculated using Vincenty’s formulae  in Microsoft Excel 2016. To test isolation by distance (IBD), pairwise FST and geographical distance (km) were examined using Mantel’s test of correlation with 10,000 permutations. Spatial autocorrelation was performed using pairwise Provesti’s genetic distance among mosquito individuals and geographic distance (km) among households with 10,000 permutations and bootstrap replications. Pairwise Provesti’s genetic distance among mosquito individuals was calculated using the ‘poppr’ package  in R . Results of the permutation were considered significant at the 5% level. In this analysis, a correlogram was produced with 45 distance classes at 1 km intervals. Both analyses yielded a correlation coefficient of the two data matrices ranging from –1 to +1, with a test for a significant relationship by random permutation. All analyses were performed using GenAlEx version 6.3 .
We observed a total of 113 alleles across 11 microsatellite loci in 21 mosquito populations from Metropolitan Manila (S3 Table). The number of alleles per loci ranged from 3 (F06) to 19 (B07) with an average of 10.25 alleles per loci, suggesting the high polymorphic nature of the selected microsatellites. A total of 91 out of 231 (39%) deviated significantly from the Hardy-Weinberg Equilibrium expectations after sequential Bonferroni correction. Seventy-two of these significant deviations indicated He > Ho, suggesting heterozygosity deficiency (S4 Table). Null alleles were present in four loci (M313, AC4, AG7 and H08) and the null allele frequency ranges from 0.00–0.33 in all loci (S5 Table). The LD test showed a total of 131 out of 1155 (11.30%) pairs of loci with significant LD after Bonferroni corrections, with no locus pair consistently correlated across all mosquito populations (S6 Table).
Table 1 shows a summary of the genetic diversity of each mosquito population. The mean number of different alleles ranged from 3.82 (TGG) to 6.36 (QZC-3) while the mean number of effective alleles areas ranged from 2.74 (QZC-1) to 3.51 (LSP-2). On the other hand, the mean allelic richness ranged from 3.24 (QZC-1) to 3.85 (MRK) while the proportion of private alleles ranged from 0.02 (LSP-2) to 0.23 (PSY). All mosquito populations except for one (MKT) did not conform to Hardy-Weinberg equilibrium expectations (He > Ho), indicating heterozygosity deficits and the possibility of inbreeding within each mosquito population.
The global FST was estimated to be 0.013 where pairwise corrected FST values among mosquito populations ranged from -0.009 to 0.064 (S7 Table). Among these comparisons, 87 out of 201 (41.4%) uncorrected pairwise FST values showed significant genetic differences (S7 Table). The effective number of migrants among mosquito populations ranged from 3.67 to 428.57 migrants (S8 Table) with 13 pairwise estimates were not calculated due to the negative or zero FST values. The results from the AMOVA showed that the FIS and FIT values were 0.131 and 0.137, respectively (both P < 0.001) (S9 Table), indicating that genetic differences exist among and within individuals.
The dendrogram based on Cavalli Sforza and Edwards distance revealed the spatial pattern and distribution of genetically similar mosquito populations (Fig 2). The “cindex” index in NbClust determined the optimal number of cluster groups to be four. One group (green line/circles, Fig 2) had 14 mosquito populations that were highly connected or genetically similar. An extended pattern can be observed from the north down to the southern parts of Metropolitan Manila.
The Mantel test between the pairwise genetic (FST) and geographic distances of all mosquito populations showed very low and non-significant correlation (R2 = 0.006, p = 0.25), indicating no isolation by distance (S1 Fig). High genetic similarity among adult mosquito individuals was only found up to 1 km based on the spatial autocorrelation analysis (S2 Fig), suggesting limited dispersal capability for Ae. aegypti.
According to Evanno et al.’s method (2005) , the most probable number of genetically differentiated groups using all individuals across mosquito populations is K = 4 (Fig 3, S3 Fig). Furthermore, the standardized datasets showed the most probably numbers of K ranged from 3 to 6 (Fig 3, S3 Fig). All barplots demonstrated no clear genetically-distinct clusters, suggesting high admixture among all mosquito populations.
Our study revealed low genetic differentiation (FST) among mosquito populations across Metropolitan Manila which correlates well with the findings from other population genetic studies of Ae. aegypti in a microgeographic-scale setting [19,21,22]. In a study conducted in Brazil, the level of genetic differentiation ranged from 0.002 to 0.094 with a maximum distance among collection sites of 30 km . Urban cities in Southeast Asian countries with spatial scales of 5–50 km also showed low genetic differentiation (0.026–0.032) . The same outcome (average FST = 0.037) was observed among villages in Thailand with geographical distances up to 27 km . Our findings, along with those of previous studies, suggest that mosquito populations within fine-scale areas consist of similar allele frequencies, exhibiting continuous and active exchange or sharing of alleles. Moreover, our study revealed the lack of isolation-by-distance (IBD) which further implies that mosquito populations are not in equilibrium of migration and genetic drift at this spatial scale.
The lack of IBD and low FST imply high gene flow and weak genetic drift among mosquito populations in Metropolitan Manila. This high gene flow may indicate the possibility of “passive” dispersal of Ae. aegypti throughout the region. Other studies claim that this mosquito vector occasionally travels long distances by taking advantage of human-aided transportation routes via land, sea, or air [26,29,62]. This is evidenced by the eggs, larvae, and adults found in commercial trucks and ships through tire importation [25,63] while larvae and pupae are littered across transportation zones such as airports  and docks/ports [29,63]. Such a mechanism could facilitate high gene flow, resulting in genetically-similar mosquito populations. Additionally, spatial autocorrelation revealed limited dispersal distance up to 1 km, suggesting the individual flight range or “active dispersal” of Ae. aegypti . The duality of both “passive” and “active” dispersal has been a common finding among human-adapted mosquitoes, especially in Ae. aegypti found in Southeast Asia [19,21,27, 65] and is considered a key factor in the persistence and resurgence of mosquito-borne diseases .
Furthermore, the low genetic differences due to weak genetic drift among mosquito populations are also reported. As evidenced by studies in Vietnam  and in the Philippines , genetic differences among mosquito populations tend to be lower during the rainy season compared to the dry season. These studies surmise that mosquito populations from the dry season recolonize breeding sites, leading to the expansion and proliferation of the mosquito vector. The large population size at this period can decrease genetic drift which leads to low genetic differences among mosquito populations. A similar phenomenon was observed in areas where insecticide is applied periodically. Mosquito populations in untreated areas facilitated high mosquito abundance and weakened genetic drift, resulting to decreased genetic differentiation among mosquito populations .
Our Bayesian analysis inferred sympatric multiple genetic clusters (K = 3 to K = 6) similar to the results reported by other population genetic studies with micro-geographic scales [19,22,67,68] and those conducted in the Philippines [28,29]. A possible explanation for this phenomenon could be the result of divergence from a single ancestry that produced multiple genetic clusters over time. It is argued that a single house or closely situated houses could act as an assembling unit in the formation of multiple genetic clusters in fine-scale areas [19,26,69]. Another possible explanation for multiple genetic clusters could be various invasions or introductions of Ae. aegypti populations from neighboring regions or provinces surrounding Metropolitan Manila. These has been shown in studies in China  and the US . We attempted to test this hypothesis by obtaining publicly available genetic data from previous studies in the Philippines [17, 28, 29], however only two microsatellite loci (AC2 and AC4) matched with our set of markers. Furthermore, we expanded our search to neighboring Southeast Asian countries (e.g. Thailand, Vietnam and Indonesia) using the database vectorbase.org , but obtained the same outcome of matched loci (AC2 and AC4). This perspective is an interesting avenue for future investigations regarding population genetic studies of Ae. aegypti in the Philippines.
Although this analysis indicated multiple genetic clusters, mixed ancestry is likewise noticeable where each individual mosquito is assigned to different genetic clusters. The bar plots showed clear patterns of genetic admixture across mosquito populations. The interpretations of our Bayesian analysis are consistent with other population genetic studies [19,22,28,29,67,68] which demonstrated low genetic differentiation with no discernable genetic structures. The constant migration of mosquito populations facilitated by passive dispersal (e.g. human-mediated transportation) could have produced the genetic admixture mosquito populations seen in our results.
The constructed dendrogram allows us to further examine the passive dispersal of Ae. aegypti by describing potential dispersal patterns and its impact on the fine-scale genetic structuring. The first group (Fig 2, green line/circles) demonstrates the spanning network of genetically similar mosquito populations from the northern up to the southern parts of the region. The long-distance genetic similarity among these mosquito populations could be due to the major roads that connect the north and south. This supports the claim that long-distance dispersal of the mosquito vector could be facilitated through human-aided transportation.
Interestingly, the mosquito populations in the eastern (Fig 2, red line/circles) and northwestern (Fig 2, yellow line/circles) areas highlight the patterns of dispersal and features that could affect fine-scale genetic structuring. These areas have a unique landscape feature of being marginally separated from neighboring cities. For example, a large highway and river split a sizeable land area in cities of Marikina (MRK) and Pasig (PSG) while large highways (e.g., expressways) divide Valenzuela (VAL) and South Caloocan (CAL-S). As a result, only a few roads and bridges connect to these cities which limits access to neighboring areas. This circumstance may limit migration between mosquito populations, resulting in formation of distinct genetic groups as shown in the dendrogram (MRK & PSG and VAL & CAL-S). Landscape fragmentation can limit gene flow among mosquito populations as evidenced in the West Indies, producing defined mosquito populations on the east and west side of a major highway . This could also explain the TGG mosquito population from Taguig City since the land area is separated by two large highways and a river which limits access to neighboring cities (Fig 2, black line/circle).
The genetic grouping shown in the dendrogram can also be explained by the occurrence and magnitude of flooding in selected cities within the region. From 2009 to 2014, Metropolitan Manila had experienced several flooding events which displaced residents and damaged properties and businesses . Porio  identified the cities of Taguig, Pasig, Marikina, Valenzuela, and Caloocan as flood-prone areas, among many others. Because of these cities marginally separated, we hypothesize that, after a flooding event, alleles within these populations could be generally preserved and its allelic combination may become distinct from other mosquito populations over time due to limited gene flow.
Our findings along with those from previous reports in the Philippines [28,29] provide a partial picture of the current population genetic status of the mosquito vector, Ae. aegypti, in the country. The low genetic differentiation observed in all the reports including ours, indicates its widespread distribution not only within cities and villages but also in different islands of the country. Human-mediated transportation (marine or road transportation) facilitate the long dispersal capability of Ae. aegypti, enabling the observed admixture within the mosquito population. This comprehensive illustration of the population biology of Ae. aegypti in the Philippines strengthens our understanding for developing better vector control strategies. If and when the country pursues innovative vector control approaches, these findings will be integral to the success of the program.