Human-Mediated Marine Dispersal Influences the Population Structure of Aedes aegypti in the Philippine Archipelago.

BACKGROUND: Dengue virus (DENV) is an extraordinary health burden on global scale, but still lacks effective vaccine. The Philippines is endemic for dengue fever, but massive employment of insecticides favored the development of resistance mutations in its major vector, Aedes aegypti. Alternative vector control strategies consist in releasing artificially modified mosquitos in the wild, but knowledge on their dispersal ability is necessary for a successful implementation. Despite being documented that Ae. aegypti can be passively transported for long distances, no study to date has been aimed at understanding whether human marine transportation can substantially shape the migration patterns of this mosquito. With thousands of islands connected by a dense network of ships, the Philippines is an ideal environment to fill this knowledge gap. METHODOLOGY/PRINCIPAL
FINDINGS: Larvae of Ae. aegypti from 15 seaports in seven major islands of central-western Philippines were collected and genotyped at seven microsatellite loci. Low genetic structure and considerable gene flow was found in the area. Univariate and multivariate regression analyses suggested that anthropic factors (specifically the amount of processed cargo and human population density) can explain the observed population structure, while geographical distance was not correlated. Interestingly, cargo shipments seem to be more efficient than passenger ships in transporting Ae. aegypti. Bayesian clustering confirmed that Ae. aegypti from busy ports are more genetically similar, while populations from idle ports are relatively structured, regardless of the geographical distance that separates them.
CONCLUSIONS/SIGNIFICANCE: The results confirmed the pivotal role of marine human-mediated long-range dispersal in determining the population structure of Ae. aegypti. Hopefully corroborated by further research, the present findings could assist the design of more effective vector control strategies.

Introduction

Aedes aegypti is the major vector of dengue virus (DENV). An estimated 96 million dengue cases and 300 million asymptomatic infections were reported globally in 2010 [1]. Currently, there is no effective dengue vaccine available and the principal way to prevent the spread of the virus and the occurrence of outbreaks is to suppress the populations of vector mosquitos. Insecticide spraying has been the most common and effective method in the past decades, but massive and indiscriminate use generated a selective pressure that led to the development of insecticide-resistant strains of mosquitos. Alternative population control methods consist in the release of modified mosquitos which will interbreed with wild populations and suppress or replace it. Depending on the technique, the released individuals are sterilized, genetically engineered or carry symbiotic bacteria of the genus Wolbachia [2].

The dispersal ability of the vector should be taken into account when implementing any of the above-mentioned population suppression techniques. If a vector can travel long distances it can quickly and effectively spread insecticide resistance mutations, symbiotic bacteria or deleterious mutations.

The flight range of Ae. aegypti is narrow (usually less than 200m, rarely up to 500m) [3-5] and its active dispersal in the wild should be limited. However, this mosquito is believed to occasionally travel long distances by taking advantage of human (land, sea or air) vehicles [6-8]. There are reports of Ae. aegypti collected from ships right after docking [9] and the trade of used tires is suspected to be one of the main contributions to the passive dispersal of Ae. aegypti immatures, both by land and sea [10]. Recently, larvae and pupae of Ae. aegypti have been for the first time collected at Narita Airport (Tokyo, Japan). It is strongly suspected that they were transported via airplane from a foreign tropical country [11]. If Ae. aegypti's migration flows are mainly driven by active dispersal, the genetic differentiation between populations should increase with the geographical distance; this is called Isolation By Distance (IBD) [12]. Otherwise, if passive human-mediated dispersal is pivotal, IBD should be weak or absent.

A remarkable amount of population genetics studies have been conducted on Ae. aegypti in many regions of the world. The genetic structure of this important vector mosquito has been analyzed at micro and macro geographical scales with different methods and purposes. Long-distance migration has been hypothesized and frequently related to human movements, but this association was not addressed in a quantitative way and only briefly brought up in the final discussion as a possible explanation of the results [13-15]. A few studies had their main focus on the effects of human transportation on Ae. aegypti's dispersal; yet, only the land transportation was taken into account [16-20]. One study undertaken in French Polynesia referred to the effect of human sea transportation; however, its main purpose was to determine if the circulation of DENV strains between Polynesian islands could be related to Ae. aegypti's migration [21]. Therefore, no study to date has been specifically aimed at quantitatively clarifying whether human-mediated long-distance marine dispersal is able to shape the genetic structure of Ae. aegypti at a macro-geographical scale.

The Philippines is an archipelago located in Southeast Asia and is endemic for dengue fever. According to the World Health Organization, the annual reported cases of dengue fever increased around 17-fold in the decade 2000–2010 [22]. Unlike other southeast Asian countries (mainly Thailand and Vietnam), very few ecological studies and no population genetics study have been conducted on Ae. aegypti in the Philippines [23-26]. A recent survey conducted by our group in the northern part of the country showed a widespread resistance to pyrethroid insecticides and a high prevalence of kdr mutations [27,28]. Furthermore, the nearly 7000 islands that make up the Philippine archipelago are connected by a dense network of cargo and passenger ships. The islands are remarkably different in size, demography and economic development; therefore, the seaports as well are strikingly different in size and commercial importance. Busy ports are more interconnected than idle ports, regardless of their geographical position. This offers a proper environment for the randomization of collection sites.

For these reasons we decided to undertake a population genetics study aimed at understanding if and how human marine transportation affects the movements of Ae. aegypti among different islands of the Philippines. Specific objectives were: (1) to analyze the genetic structure of Ae. aegypti among different islands; (2) to quantitatively determine if natural (geographical distance) or anthropic (port size and intensity of marine transportation) factors affect the observed genetic structure.

Methods

Sampling

In order to analyze the effect of human sea transportation on the migration of Ae. aegypti, seaport locations on several islands were sampled, as populations of mosquitos residing near the seaports are more likely to end up resting or laying eggs on the ships (either cargo or passenger). Each site corresponded to one seaport and all individuals from a site were considered as one separate population, therefore the terms "site", "population" and "port" will be used interchangeably, hereinafter. The survey was conducted between September and October 2013 in the central-western area of the Philippines (Region IV-A and IV-B). Larvae were collected from 15 ports of different levels of connectivity and development, across seven islands of different sizes and degrees of urbanization (Fig 1 and Table 1). At each site, as many water-filled containers as possible (artificial or natural) were inspected in a 1.5 km-radius area; the radius was calculated from the center of the port. In case the container was located in a private area, the inspection was conducted only with consent of the owner. Collected larvae were conserved in absolute ethanol and identified by microscopic analysis [29,30] or by PCR amplification with species-specific primers [31]. Only larvae identified as Ae. aegypti were included in the genetic study. In order to reduce the proportion of sibling individuals (which overestimate the genetic difference between populations), Ae. aegypti's flight range and its so-called "skip oviposition" behavior were taken into account [32] and a minimum distance of 200m was kept between the containers; then, one out of five larvae were randomly selected from each container and included in the genetic analysis. For port 6 one out of six larvae were selected. In case it was not possible to keep the minimum distance of 200m, larvae from multiple containers were grouped as if collected from the same.

Fig 1

Map of the sites where Ae. aegypti were collected.

For detailed populations information refer to Table 1.

Table 1

Details of the populations studied.

Site Number 1	Site Name	Latitude	Longitude	Containers 2	Individuals 3	N 4
1	Ambulong	12°29'34.10"N	122°29'27.78"E	8/10	346	67
2	Balanacan	13°32'0.02"N	121°51'59.28"E	7/12	112	24
3	Batangas	13°45'5.68"N	121°2'36.99"E	8/9	288	56
4	Brooke's Point	8°46'16.83"N	117°50'9.97"E	4/11	235	49
5	Cajidiocan	12°22'2.06"N	122°40'41.43"E	6/11	268	53
6	Calapan	13°25'37.12"N	121°11'46.03"E	7/11	301	48
7	Cawit	13°22'55.07"N	121°49'27.23"E	6/13	60	14
8	Lucena	13°54'18.18"N	121°37'18.06"E	7/12	279	56
9	Odiongan	12°25'4.53"N	121°59'46.95"E	4/12	195	40
10	Puerto Princesa	9°44'22.93"N	118°44'30.24"E	8/14	204	43
11	Romblon	12°34'52.48"N	122°15'59.25"E	5/15	93	19
12	Roxas	12°35'37.56"N	121°31'10.47"E	5/13	140	30
13	San Agustin	12°34'3.81"N	122°8'5.02"E	4/6	324	64
14	San Jose	12°19'51.36"N	121°5'12.97"E	7/15	152	31
15	Torrijos	13°19'21.40"N	122°5'15.79"E	8/12	231	48

1Refer to Fig 1 for geographical location of each site.

2Containers inspected (positive for Ae. aegypti/total).

3Total Ae. aegypti collected.

4Individuals selected for the study.

Genotyping

Genomic DNA was extracted from individual larvae with RedExtract-N-AmpTM Tissue PCR Kit (SIGMA, USA) following the manufacturer's protocol. Seven microsatellite loci were scored for each individual. Each PCR reaction (10μl) contained 400nM of M13 fluorescent primer (FAM-labeled) [33], 200nM of M13-tailed forward primer, 200nM of reverse primer, 1X of Buffer (Mg2+ 2mM included), 200nM of each dNTP, 0.25U of Takara ExTaq (Takara BIO, Inc.) and 1μl of sample. Thermocycling conditions were: (5' at 94°C) for one cycle, (30'' at 94°C, 30'' at Ta, 40'' at 72°C) for 35 cycles, (5' at 72°C) for 1 cycle and 4°C as holding T. T of annealing (Ta) was specifically determined for each pair of primers (S1 Table). After PCR amplification, 1.5μl of product was added to 8.5μl of a GeneScan™ 500 ROX Size Standard:formamide solution (1:100), heated at 95°C for 5' and run on an Applied Biosystems 3730xl DNA Genetic Analyzer. The allele scoring was conducted with the GeneMapper Software 5 (Applied Biosystems).

Marker selection

Seven microsatellite markers were selected out of 74 loci previously characterized in other studies [34-38]. After a search in the database, 33 were found to reside on different supercontigs. Preliminary PCR reactions with different Ta were conducted in order to find the most suitable for each primer pair (S1 Table). Trial runs on capillary electrophoresis were performed to select the markers that reliably amplified in most of the samples, with sufficient allelic variability and that could be clearly and easily scored in the software GeneMapper. Many loci were discarded because the size of some of their allele peaks were misaligned with the sizes expected from their repetitive unit size; most of those loci turned out to be located near poly-A segments or to have deletions/insertions not related to the short tandem repeat of interest.

Data analyses

The software GENEPOP v4.2.1 [39,40] was used to compute the basic population genetic parameters. The Inbreeding Coefficients (Fis) were calculated for all loci across populations according to Weir and Cockerham [41], in order to detect any statistically significant deviation from Hardy-Weinberg Equilibrium (HWE). Each pair of loci was also tested for Linkage Disequilibrium (LD). Markow Chain parameters were set at 10,000 dememorizations, 20 batches, 5000 iterations for the former analysis and at 10,000 dememorizations, 100 batches, 5000 iterations for the latter one. GENEPOP was also used to calculate total and pairwise Fixation Indexes (Fst); P-values for pairwise Fst were computed with the software Arlequin v3.5.1.3 [42].

To test the hypothesis that human transportation influences mosquitos' dispersal, eight variables were elaborated; two were obtained from satellite images and six from official data published online by the Philippine Statistics Authority (Table 2 and S2 Table). The variable named "Distance" accounts for the geographical location of the ports. The other seven variables account for port size and degree of ship connectivity, therefore we will refer to them as Anthropic Variables (AVs) hereinafter. Two AVs ("Inhabitant" and "Density") are demographic data on the municipalities where ports are located. The variable named "Dock" expresses the physical size of a port, by measuring the perimeters of its concrete docks. The last four AVs ("Vessel", "Tonnage", "Cargo", "Passenger") measure the intensity of ship traffic at ports. The statistics for the last four AVs were incomplete and their correspondent values for ports 1, 5, 13 and 15 were missing. In order to interpolate them, four bivariate linear regressions with Inhabitant and Dock as independent variables were run. Inhabitant and Dock were chosen because of high correlation with Vessel, Tonnage, Cargo and Passenger (Pearson's r≥0.47, S3 Table). The command "lm" in the software R for MachIntosh (v3.1.1 GUI 1.65 Snow Leopard build) was used for the purpose.

Table 2

Definitions of the variables used to test the hypothesis that human transportation influences population structure of Ae. aegypti in the central-western Philippines.

Name	Definition
Fst	Fixation Index calculated for each pair of Ae. aegypti populations, using allelic frequencies at seven microsatellite loci (S1 Table); high Fst values indicate genetic differentiation, while low Fst indicate genetic similarity
Distance	Euclidean distance in km between each pair of ports 2
Inhabitant 1	Total human population of the municipality to which each port belongs 3
Density 1	Total human population of the municipality to which each port belongs, divided by the total area of the municipality in square kilometres 3
Dock 1	Total sum of the perimeters (in meters) of all the concrete docks of a port 2
Vessel 1	Total annual number of domestic vessels visiting a port, inbound and outbound (average from 2000 to 2011) 4
Tonnage 1	Total annual Gross Register Tonnage (a ship's total internal volume expressed in register tons) of the domestic vessels visiting a port, inbound and outbound (average from 2000 to 2011) 4
Cargo 1	Total annual domestic cargo throughput of a port, in metric tons, inbound and outbound (average from 2000 to 2011) 4
Passenger 1	Total annual number of passengers visiting a port, inbound and outbound (average from 2000 to 2011) 4

1Originally single data for a single port, then transformed to pairwise using the formula a+b where a and b are single values for each port.

2Calculated from satellite images using the "ruler" function of Google Earth.

3Official data from the Philippine Statistics Authority—National Statistics Office (http://web0.psa.gov.ph/).

4Official data from the Philippine Ports Authority (http://www.pdosoluz.com.ph/).

The command "princomp" in R was used to run a Principal Component Analysis (PCA) on the AVs. As the data showed considerable heteroscedasticity they were log-transformed before PCA [43].

As the AVs were expressed as a single value for each port (15 values per variable), while Fst values were pairwise (105 total port-by-port pairs, excluding same-port pairs), in order to quantitatively compare them, the values for the AVs were transformed to pairwise; the formula a+b was used, were a and b are the corresponding value for each two ports. Pearson correlation coefficients were computed in R with the command "rcorr". To detect collinearity, the variance inflation factors (VIF) between eight predictors were computed (command "vif"). The command "lm" was used to run the multivariate linear regression models. The best model was chosen with the command "step", which performs a stepwise exclusion of each predictor and calculates the Akaike Information Criterion. Simple and partial Mantel tests were run using the online software "Isolation by distance web service" (IBDWS) [44]. One-tailed P-values were calculated with 10,000 bootstrap randomizations.

The Bayesian clustering method implemented in the software STRUCTURE v. 2.3.4 [45] was used to assign individual Ae. aegypti to genetic clusters based on the microsatellite allelic frequencies. Admixture was selected as ancestry model with LOCPRIOR function. Alpha was allowed to vary and was separately calculated for each population. A correlated allele frequency model was chosen, with lambda set to one. Ten independent replicates were run with a burn-in step of 300,000 iterations followed by 500,000 replications. Average L(K), ΔK [46] and proportions of membership (i.e. the probability of a population to belong to a certain cluster, we will hereinafter refer to it as POM) were computed in Excel. Pearson correlation coefficients between POM and principal component 1 (PC1) were computed in R with the command "cor.test".

Results

An average of 211 Ae. aegypti larvae per site were collected (total = 3166). At each site between six and 15 containers were inspected and between four and eight were positive. One out of five larvae was randomly selected from each container and an average of 43 per population was genotyped (total = 642) (Table 1).

Marker validation

The test for deviation from HWE calculated 105 Fis values [41], one for each locus-population pair (S4 Table). Nine (8.6%) showed significant departure from HWE (after Bonferroni sequential correction), eight towards a heterozygosity deficit and one towards a heterozygosity excess. However, no locus showed significant values at more than three populations.

The test for LD calculated a total of 315 values (21 locus-by-locus pairs for each population), of which 28 (8.9%) remained significant after Bonferroni sequential correction. However, no locus-by-locus pair showed significance consistently across populations. Therefore this significance was unlikely due to linkage. None of the markers was excluded from the analysis.

Genetic structure of Ae. aegypti

For most loci between four and six alleles were detected; only AG5 had 11. The Fst for all loci and populations was 0.056, while the pairwise Fst ranged from 0.005 and 0.147; after Bonferroni correction 79/105 (75.2%) remained significant (P<0.0005) (Table 3).

Table 3

Fixation indexes (Fst) calculated for each pair of populations studied.

	1	2	3	4	5	6	7	8	9	10	11	12	13	14
2	0.060*
3	0.051*	0.010
4	0.119*	0.029*	0.069*
5	0.099*	0.045*	0.049*	0.087*
6	0.081*	0.018	0.023*	0.077*	0.040*
7	0.077*	0.019	0.062*	0.083*	0.069*	0.060*
8	0.056*	0.023	0.009	0.075*	0.041*	0.038*	0.073*
9	0.147*	0.042*	0.051*	0.072*	0.105*	0.066*	0.131*	0.066*
10	0.063*	0.021	0.016	0.057*	0.046*	0.059*	0.074*	0.022*	0.077*
11	0.024	0.019	0.022	0.073*	0.070*	0.061*	0.043	0.022	0.081*	0.029
12	0.047*	0.007	0.009	0.047*	0.044*	0.037*	0.034	0.017	0.079*	0.023	0.005
13	0.103*	0.059*	0.049*	0.075*	0.118*	0.117*	0.086*	0.067*	0.118*	0.062*	0.051*	0.025
14	0.080*	0.021	0.019	0.097*	0.031*	0.020	0.053*	0.029*	0.067*	0.029*	0.030	0.034	0.104*
15	0.106*	0.039*	0.034*	0.093*	0.064*	0.061*	0.057*	0.051*	0.101*	0.035*	0.049*	0.046*	0.081*	0.023

*P<0.0005.

Statistical analysis of factors affecting Ae. aegypti's genetic structure

Pearson correlation coefficient for Distance against Fst was close to zero and non-significant, while it was highly significant for all the AVs against Fst (P<0.01) (Table 4). Tonnage and Passenger were excluded from the multivariate linear regression because of high collinearity (square root of VIF>10). Distance, Inhabitant, Dock and Vessel were dropped after stepwise analysis. The multiple R2 was 0.24 (P<0.001). The regression coefficients for Density and Cargo indicated a negative association with pairwise Fst (P<0.001) (Table 4 and Fig 2). The square roots of VIFs are displayed in S5 Table.

Table 4

Results of regression models.

Predictor	Pearson's r	Multiple linear regression
		Excluded 1	Dropped 2	Significant 3
Distance	0.09		✗
Inhabitant	-0.44		✗
Density	-0.34			-1.03e-05
Dock	-0.36		✗
Vessel	-0.29		✗
Tonnage	-0.3	✗
Cargo	-0.37			-4.16e-08
Passenger	-0.33	✗

In all cases pairwise Fst is the dependent variable. Bold indicates P<0.01 for Pearson and P<0.001 for multiple linear regression.

1Predictors excluded before running the model because the square root of Variance Inflation Factor (VIF) was >10.

2Predictors included in the model, but dropped by the stepwise analysis.

3Regression coefficients of the predictors confirmed by the model. Refer to Table 2 for detailed description of the predictors.

Fig 2

Three-dimensional scatter plot with regression plane that visualizes the effect of the variables called “Density” and “Cargo” on pairwise Fst values.

Both predictors are negatively correlated to Fst and the relative regression coefficients are very low in absolute value. This means that a remarkable increase in Cargo and Density is necessary to cause a unit decrease in Fst values. The interpretation is that where cargo shipments are intense and human population is dense Ae. aegypti are more genetically similar. This suggests an influence of human transportation on mosquitos’ migration. Refer to Table 2 for detailed description of the variables.

The PCA on the seven AVs returned a PC1 able to explain 74% of the total variance of the dataset; the directions of the vectors representing the AVs were almost opposite to the PC1 (Fig 3). It is therefore reasonable to assume that higher PC1 values indicate ports located in less developed and less connected municipalities and vice versa.

Fig 3

Principal component analysis used to summarize seven variables related to port size and connectivity; plot of the first two PCs.

Red arrows represent the vectors of the seven variables and each point in the plot is a port (Fig 1 and Table 1). The direction of principal component 1 (PC1) is almost opposite to the variables and it explains 74% of the total variance in the dataset; therefore, ports with low values of PC1 were considered big and highly connected (and vice versa). Refer to Table 2 for a detailed description of each variable.

STRUCTURE analysis could detect no genetic structure unless the option LOCPRIOR was chosen. This function uses information on sampling locations to assist the clustering in case the signal of structure is too weak [47]. A considerable degree of admixture was detected and most of the individuals were assigned to several clusters with almost equal probability. Nevertheless, with increasing K some populations had increasing probability to be assigned to a specific cluster (Fig 4). Following the method proposed by Evanno et al. [46] we found peaks of ΔK at K = 3, K = 6 and K = 9 the last one being the highest (S1 Fig). The peak at K = 3 was higher than K = 6, but the bar plot of the latter was considerably more informative. The highest ΔK was at K = 9, but L(K) decreases considerably between K = 6 and K = 9 (S1 Fig) and the additional structure explained by K = 9 is not remarkable (Fig 4). Therefore, we discuss clusterings at both K = 6 and K = 9. Pearson correlation coefficients of POM regressed on PC1 scores were both high and significant (for K = 6, r = 0.76, P<0.01; for K = 9, r = 0.57, P<0.05).

Fig 4

STRUCTURE bar plots.

Each individual is represented by a vertical bar whose colors show the probability to be assigned to specific clusters. Populations are separated by vertical black lines and identified by the numbers at the bottom (Fig 1 and Table 1). K = number of genetic clusters. Overall, a low level of genetic structure was detected. At K6 and K9 most populations from busy ports show genetic admixture, while most populations from idle ports are genetically clustered (compare with Fig 3).

Discussion

The results suggested low genetic structure and considerable gene flow. Pairwise Fst values were generally low (0.005~0.147), in accordance with other studies conducted in Florida (0.004~0.171), Pakistan (0.013~0.17) and Southeast Asia (0.02~0.185) [17,19,48]. In the study area the migration of mosquitos was not affected by geographical distance, but increased with the intensity of human marine transportation. In our understanding, this is the first study that qualitatively and quantitatively shows how human ship transportation affects the migration of Ae. aegypti.

Regression models

Both univariate and multivariate regression models excluded IBD. In fact, Distance/Fst Pearson's r was almost zero (0.09; non-significant) and Distance was dropped by multiple linear regression stepwise analysis (Table 4). On the other hand, Pearsons supported a negative correlation between Fst and the AVs (all r values ≤-0.3; P<0.01). The multivariate analysis confirmed the correlation only for Density and Cargo (P<0.001). The multiple linear regression coefficients were negative and very small in absolute value, therefore pairwise Fst decreases very slowly as Density and Cargo increase (Table 4); this is also evident in the scatter plot (Fig 2).

Our interpretation is that increasing cargo shipments increase the probability for Ae. aegypti to be accidentally transported, decreasing the genetic differentiation and consequently the Fst. Similarly, more densely populated municipalities have likely busier ports, with more congested ship traffic and offer more opportunities to Ae. aegypti to be transported.

The results also suggest that cargo shipments transport Ae. aegypti more effectively than passenger ships. If we consider that cargo ships (like barges) are usually bigger in size than passenger ships, they may offer a higher number of breeding sites (i.e. collections of water). Moreover, the cargo storage compartments are usually kept closed during the navigation, while passenger ships' decks are often not enclosed; this could cause more adult mosquitos to enter and being effectively transported by the cargo ships. Another hypothesis (not excluding the previous ones) is that cargo shipment routes are implemented earlier than passenger ships' routes during the development of the ports. In other words, ports of small and intermediate sizes could be connected by more cargo than passenger ships, while busiest ports operate a similar amount of the two ship types. This hypothesis is supported by the scatter plots between the four "ship connectivity" predictors (i.e. Vessel, Tonnage, Cargo and Passenger); all of them are highly correlated, but only Cargo shows a non-linear relation with the others. Values of Cargo increase earlier than values of the other three predictors (S2 Fig).

Bayesian clustering

STRUCTURE analysis showed overall a considerable admixture and most of the genetic clusters were shared by different populations. Nevertheless, the populations showed variable degrees of genetic structure, from totally admixed (port 3) to remarkably isolated (port 1) (Fig 4). At both K = 6 and K = 9, the level of admixture of the populations was strongly correlated to the development of the corresponding port (POM regressed against PC1). In other words, mosquitos from busier ports were more genetically admixed and mosquitos from smaller ports were more genetically isolated.

Interestingly, the individuals had mixed ancestry while the populations as a whole were homogeneous (this produced a pattern of "horizontal stripes" in the bar plot) (Fig 4). This could be interpreted as the result of multiple connections between many of these ports for relatively long time. This pattern is different from other studies that analyzed Ae. aegypti's population structure at broader spatial scales (several countries and continents); in these studies the ancestry of individual mosquitos was more definite, but mosquitos from the same geographical populations were assigned to different genetic clusters (this produced "vertical stripes" in the STRUCTURE bar plot) [49,50]. A reasonable explanation is that mosquitos migrate more frequently and are more admixed within the same country than between countries and continents; as mosquitos are more distantly related at broader geographical scales, migrant individuals can be more clearly identified.

Limitations

Main limitation of this study could be the sample size (Table 1), as some of the populations were analyzed with few individuals (especially ports 7 and 11); this could bias the estimation of Fst. Most locations yielded more than 100 individuals, but from few breeding sites. Therefore, we had to compromise between reducing sibling individuals (which could also bias the analysis) and have enough specimens to genotype. After Bonferroni correction 26 pairwise Fst values (24.8% of the total) were not statistically significant and most of them corresponded to the ports with smaller sample size (Table 3). On the other hand, only four (3.8% of the total) remained non-significant if Bonferroni correction was not applied (S1 Dataset). Besides, simulation showed that between 25 and 30 individuals are enough to correctly estimate the allele frequencies with five-nine microsatellite markers [51].

The data from the Philippines Statistics Authority used to elaborate Vessel, Tonnage, Cargo and Passenger was incomplete and it was necessary to estimate the corresponding values for four ports.

However, we believe that these limitations did not undermine the general meaning of our findings.

Conclusion

Our overall interpretation is that the human ship transportation substantially shapes the dispersal patterns of Ae. aegypti in the central-western part of the Philippines. Regardless of their geographical location, large ports in highly-populated municipalities exchange individuals of Ae. aegypti at a faster rate than small-sized ports in less developed municipalities. Insecticide spraying, release of sterile males, release of Wolbachia symbionts: whichever vector control strategy will from now on be applied in the Philippines and other archipelagoes, these results should be taken into account.

L(K) and ΔK at different values of K.

Peaks of ΔK were found at K3, K6, K9. L(K) were averaged over 10 independent replicates.

(PDF)

Click here for additional data file.

Scatter plot matrix of the predictors used in univariate and multivariate regressions.

Density is non-linearly correlated to predictors related to port size (i.e. Inhabitant and Dock), while Cargo is non-linearly correlated to predictors related to ship connectivity (i.e. Vessel, Tonnage and Passenger). Then, Density and Cargo are totally non-correlated between each other. For detailed description of the predictors refer to Table 2.

(PDF)

Click here for additional data file.

List and characteristics of the microsatellite markers used in the study.

(DOCX)

Click here for additional data file.

Summary statistics of the predictors used in univariate and multivariate regression.

(DOCX)

Click here for additional data file.

Pearson correlation coefficients between the predictors used in univariate and multivariate regression.

(DOCX)

Click here for additional data file.

Inbreeding coefficients (Fis) for each microsatellite locus across each population.

(DOCX)

Click here for additional data file.

Square roots of variance inflation factors (VIF) of each predictor in two multiple linear regression models.

(DOCX)

Click here for additional data file.

Scanned copy of the result of the ethical examination.

(PDF)

Click here for additional data file.

Zipped folder containing the dataset organized in separate Excel sheets.

(ZIP)

Click here for additional data file.