The genetic structure of Aedes aegypti populations is driven by boat traffic in the Peruvian Amazon

In the Americas, as in much of the rest of the world, the dengue virus vector Aedes aegypti is found in close association with human habitations, often leading to high population densities of mosquitoes in urban settings. In the Peruvian Amazon, this vector has been expanding to rural communities over the last 10–15 years, but to date, the population genetic structure of Ae. aegypti in this region has not been characterized. To investigate the relationship between Ae. aegypti gene flow and human transportation networks, we characterized mosquito population structure using a panel of 8 microsatellite markers and linked results to various potential mechanisms for long-distance dispersal. Adult and immature Ae. aegypti (>20 individuals per site) were collected from Iquitos city and from six neighboring riverine communities, i.e., Nauta, Indiana, Mazan, Barrio Florida, Tamshiaco, and Aucayo. FST statistics indicate significant, but low to moderate differentiation for the majority of study site pairs. Population structure of Ae. aegypti is not correlated with the geographic distance between towns, suggesting that human transportation networks provide a reasonable explanation for the high levels of population mixing. Our results indicate that Ae. aegypti gene flow among sub-populations is greatest between locations with heavy boat traffic, such as Iquitos-Tamshiaco and Iquitos-Indiana-Mazan, and lowest between locations with little or no boat/road traffic between them such as Barrio Florida-Iquitos. Bayesian clustering analysis showed ancestral admixture among three genetic clusters; no single cluster was exclusive to any site. Our results are consistent with the hypothesis that human transportation networks, particularly riverways, are responsible for the geographic spread of Ae. aegypti in the Peruvian Amazon. Our findings are applicable to other regions of the world characterized by networks of urban islands connected by fluvial transport routes.

Introduction Anthropogenic activities such as trade and transportation contribute to the unintentional spread of invasive organisms across the globe [1,2], resulting in serious consequences for public health, agriculture, the economy, and native ecosystems [3][4][5]. Pathogens and their insect vectors are important examples of invasive organisms that directly impact human health.
The invasive mosquito, Aedes aegypti, is the primary vector of dengue, urban yellow fever, and Zika viruses, is an important vector of chikungunya virus [6][7][8], and is a competent or suspected vector of Mayaro virus [9]. Thought to be African in origin, Ae. aegypti most likely was transported to the Americas via ships used for the transport of slaves and goods in the 15 th -19 th centuries [10][11][12]. Ae. aegypti was first reported in Peru in 1852, but was declared eradicated in 1958 following the success of a large-scale Pan American Health Organization yellow fever control program [13]. The Amazonian city of Iquitos was the first documented site of Ae. aegypti (in 1984) and dengue (in 1990) reestablishment in Peru [14,15]. In recent years, Ae. aegypti mosquitoes have been expanding geographically from urban to peri-urban and rural areas throughout Peru and the rest of Latin America [16][17][18].
Ae. aegypti dispersal can occur in one of two ways: the slower, active dispersal of flying adult females in search of bloodmeals or oviposition sites, or the faster, passive human-mediated dispersal, by which humans unintentionally transport mosquitoes via vehicle traffic (boats, cars, planes, etc). The latter can involve transport of eggs, larvae, pupae, or adults, all of which have been documented in vehicles, particularly in boats [10,[19][20][21][22]. Our previous research in the Peruvian Amazon demonstrated Ae. aegypti infestation on different vehicles commonly used for trade and transportation, including large barges, medium-sized barges, and buses [23,24]. Characterizing the relative role and importance of both active and passive dispersal mechanisms is paramount for understanding vector population structure and the dynamics of pathogen transmission.
Previous studies have characterized Ae. aegypti population genetic structure at global scales [25][26][27][28][29], in regions within countries [30][31][32], and within cities [33,34]. Ae. aegypti genetic differentiation across various spatial scales is likely due to its limited flight range (rarely exceeding 100m under natural conditions [35][36][37]), heterogeneities in insecticide application, human population densities, and water storage habits [33,38,39]. At coarser scales (i.e., between nations, or regions within nations), human transportation networks have been implicated as a driver of Ae. aegypti gene flow; that is, relatedness between populations has been shown to correlate with major highways and waterways [25,32]. The Amazonian city of Iquitos, Peru R01 AI069341-01 and P01 AI098670 (TW Scott, PI) and by the Research and Policy for Infectious Disease Dynamics (RAPIDD) program of the Science and Technology Directory, Department of Homeland Security, and Fogarty International Center, National Institutes of Health. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
presents is a unique setting in which to address questions of Ae. aegypti gene flow and human transit, as people in this region rely heavily on rivers as a means of transport between locations. Iquitos is the site of many longitudinal studies related to dengue epidemiology and Ae. aegypti ecology [40], yet little is known about its population structure in this region. In this study, we used Ae. aegypti samples collected from several sites within and around Iquitos to characterize the population structure, and preliminarily explored whether genetic relatedness among mosquito populations is driven by human transportation networks.

Ethics statement
Permission for this study was granted by the Loreto Regional Health Department, and the study protocol was approved by the NAMRU-6 Institutional Review Board in compliance with all applicable federal regulations governing the protection of human subjects (protocol number NAMRU6.2012.0039). In addition, the Emory University Institutional Review Board determined that this study does not represent human subjects based research.

Study area
Accessible only by plane or boat, Iquitos is surrounded by other, smaller settlements that are primarily connected to one another via river networks. The only significant source of terrestrial transit is the 95km Iquitos-Nauta Highway, connecting Iquitos (population, pop henceforth: 406,340) to the smaller city of Nauta (pop: 13,983) [41]. Once the epicenter of the rubber industry in the early 1900s, the Iquitos economy now relies predominantly on oil and timber exportations in addition to tourism. This setting is ideal for studying the invasion dynamics of Ae. aegypti, because the region's inhabitants are dependent on both fluvial and terrestrial routes for trade and transportation.

Mosquito collections
Mosquitoes were collected from the city of Iquitos, and the neighboring towns of Nauta (pop: 13 [42] using shapefiles generated from prior research activities [43].) After asking permission to survey the household for Ae. aegypti mosquitoes, team members collected adult and immature mosquitoes either through aspiration of adults or larval surveys. Mosquito collection methods are described in detail elsewhere [44]. Within communities, we diversified mosquito sampling to the extent possible: we collected a maximum of 3 mosquitoes from any given household, and sampled 1 out of every 10 houses. More than 20 individuals were collected for each sampling location, although mosquitoes collected from Indiana (n = 18) and Mazan (n = 14) were grouped into one population to ensure adequate sample sizes for the calculation of inbreeding coefficient F IS , F ST , and Bayesian clustering analysis (described in detail below). Indiana and Mazan are only 1.2 km apart and our data indicated that mosquito populations in those two locations were genetically indistinguishable (F ST = 0.032). The Euclidean distances between pairs of towns ranged from approximately 15km (Aucayo and Tamshiaco) to 125km (Nauta-Indiana/Mazan). (S1 Table  demonstrates Euclidean pairwise distances between towns. ) We collected mosquitoes from multiple locations in Iquitos to determine whether mosquitoes in port areas were more closely related to those in the surrounding towns (Fig 1). The location Iquitos A (Puerto Masusa) serves as a major hub of fluvial transit and predominantly  harbors large barges, which carry cargo and passengers throughout the Peruvian Amazon distances up to~500km. Iquitos B (Huequito) is a secondary fluvial port, primarily harboring medium-sized barges that also carry cargo and passengers, but that travel more locally (up to distances of~250km) [24]. We also sampled mosquitoes from the interior of Iquitos (Iquitos C, district of Iquitos) and from a recently urbanized southern neighborhood near the start of the Iquitos-Nauta highway (Iquitos D, district of San Juan).

DNA extraction and amplification
Methods used for genotyping mosquitoes with microsatellite markers are described by Wong et al [45]. Genomic DNA from the mosquito body was purified by potassium acetate/ethanol precipitation [46]. In a multiplex polymerase chain reaction (PCR), we amplified 8 previously described microsatellite markers [47,48]. PCR products were diluted 1:60 or 1:40 in ddH 2 O, and submitted to the University of California Davis College of Agriculture and Environmental Sciences Genomics Facility for fragment analysis on an ABI 3730 XL capillary electrophoresis sequencer (Life Technology Corp.). GS600 LIZ size standard (Life Technology Corp.) was included with each sample to determine the size of individual peaks. ABI Peak Scanner software (Applera Corp., Norwalk, CT) was used to visualize resulting chromatograms. After identifying fragments, alleles were assigned using the MsatAllele package in R [49,50].

Microsatellite data analysis
We used the microsatellite data to calculate observed and expected heterozygosity values, in addition to inbreeding coefficients, F IS , for each site and genetic locus in Arlequin v3.5.2.2 [51]. We also tested for departure from Hardy-Weinberg equilibrium (exact test with 1,000,000 Markov chain steps and 100,000 dememorization steps) in Arlequin v3.5.2.2. A Sidak correction for multiple comparisons (for a total of 72 tests) was applied to determine significant deviation from Hardy-Weinberg equilibrium (p < 0.00071). F ST was calculated in Arlequin v3.5.2.2 (10, 000 permutations, significance level p < 0.05) [51]. We first conducted a pairwise F ST analysis on all nine sampling locations to understand the degree of genetic differentiation among all sample locations including Iquitos. We then conducted a second pairwise F ST analysis combining Iquitos samples to understand degree of connectivity between the city center and the five villages outside the city. Mantel tests conducted in the package ade4 in R [50,52] were used to test genetic isolation (10,000 permutations) by three different distance models for the data that combined Iquitos samples: 1) Euclidean distance, the shortest straight-line distance between two locations, 2) fluvial path distance, the river route between towns, and 3) shortest path distance, the shortest accessible fluvial or terrestrial route between two locations. Input F ST values used to generate isolation by distance plots are shown in S2 Table. We also developed a "Propagule Pressure Index," combining the probability of Ae. aegypti infestation in different vehicle types with the frequency of travel between Iquitos and surrounding towns. We calculated infestation probabilities from data collected in 2013 across six different vehicle types common in Iquitos including: large and medium size barges, water taxis, speed boats, buses, and taxis [24]. Our results indicated that some vehicle types were consistently infested with Ae. aegypti across multiple months (71% of large barges, 35% of medium-sized barges, and 12.5% of buses). Simultaneous with entomological surveys, we interviewed vehicle drivers to determine the frequency of travel between Iquitos and surrounding towns for each vehicle type to estimate the approximate number of vehicles traveling to each town. The Propagule Pressure Index (PrPI) is calculated as follows: Where: γ i = probability of vehicle infestation θ i = number of trips from Iquitos to surrounding towns i represents individual vehicles of a certain type j represents different vehicle types (large barges, medium barges, and buses) PrPI therefore represents the number of trips made between pairs of locations, weighted by the probability of Ae. aegypti infestation for each vehicle type. Values of PrPI were plotted against F ST to visually explore the relationship between these two variables. However, no statistical tests were conducted due to low sample sizes-pairwise calculations of PrPI were only available between Iquitos and other towns, leaving a total of five observations.
We also estimated the genetic divergence of the populations by an analysis of molecular variance (AMOVA) in Alrequin v3.5.2.2 [51]. The geographic structure considered in this analysis is shown in Table 1; Aucayo, Barrio Florida, Indiana-Mazan, and Tamshiaco were counted as separated populations, but Iquitos sites were given hierarchical consideration. The total variation observed was attributed to differences between individuals within populations, and among populations. Accordingly, we calculated pairwise F statistic analogues characterizing the variation among individuals and populations (10,000 permutations, significance level p < 0.05).
We investigated the regional population structure of Ae. aegypti with the program Structure v. 2.3.4 [53]. Structure uses estimated allele frequencies to compute the likelihood that a given genotype originated from a genetic cluster. The result is that probabilistic estimates of population membership coefficients are assigned to each individual. Ancestral genetic admixture within an individual is observed when an individual has more than one population group is assigned. (For example, individual A is genetically admixed if the probability of belonging to group 1 is 0.6 and to group 2 is 0.4-the sum of these components is 0.6 + 0.4 = 1.) In Structure, we used an admixture model with uncorrelated allele frequencies to avoid the risk of overestimating the number of populations (100,000 burn-ins and 200,000 Markov Chain Monte Carlo runs after the burn-in period). We started simulations with K = 10, to allow for the possibility of more genetic clusters than sampling locations (9 sampling locations including 4 Iquitos locations, and 5 towns), and then ran simulations for K values of 10 through 1. For each K, we ran 10 simulations to ensure consistency between runs, and used the log likelihood [53] and DeltaK method [54] to determine the most likely number of genetic clusters.

Results
After applying a Sidak correction for multiple comparisons (p < 0.00071) for all populations, significant heterozygous deficits were detected in eight comparisons (both statistically significant and F IS > 0.50) among four loci (B07, AC1, AG2, H08) ( Table 2). These comparison were distributed across five sites including Barrio Florida (B07), Tamshiaco (B07), Iquitos A (AC1 and AG2), Iquitos C (B07 and H08), and Iquitos D (AC1 and B07). At a site-level, significant and low to moderate levels of heterozygous deficits were observed in Nauta, Tamshiaco, and all sites within Iquitos (F IS range among these sites: 0.19629 to 0.37096).Of 72 tests, 45 (62.5%) loci from all populations were found to be in Hardy-Weinberg equilibrium. The majority of deviations were from Nauta and Tamshiaco; when those populations were excluded, 73.2% of loci (41 of 56) were in Hardy-Weinberg equilibrium.
Pairwise F ST values demonstrated low to moderate differentiation for the majority of site pairs ( Table 3). Mosquitoes from Barrio Florida and Nauta were significantly differentiated Our analysis of isolation by distance showed no clear relationship between genetic distance (F ST ) and geographic distance for the three standard distance models (Fig 2). We were unable to statistically evaluate the relationship between PrPI and genetic distances because of insufficient transportation network data. AMOVA results showed that 71.5% of the total variation observed was attributable to differences within individuals; only 4.1% of the total variation was attributable to differences among samples collected within Iquitos, indicating a lack of genetic structure within the city.
Results from the DeltaK analysis of Structure output indicated that there were most likely three genetic clusters (S1 Fig). That is, the DeltaK statistic was highest for k = 3 groups (Del-taK = 9.5) in comparison with all other possible numbers of groups (DeltaK range: 0.18 to 1.6 for k = 1, k = 2, and k = 4 through 10). Structure results showed a clear pattern of genetic admixture from all the sample populations, suggesting steady gene flow within this region (Fig  3). Of 339 individual mosquitoes, 181 (53.4%) showed dominant haplotype membership associated with a single genetic group (represented in blue in Fig 3). Dominant haplotype membership (> 0.8) to the same group (depicted in blue in Fig 3) was observed among 84.4% of mosquitoes collected from Indiana/Mazan, 87.5% of mosquitos from Iquitos B, 75% of mosquitoes from Aucayo, and 64.5% of individuals from Barrio Florida. Mosquitoes collected from Tamshiaco and Nauta showed the greatest degree of ancestral admixture, with each of the three genetic groups represented approximately equally.  [EXSCINDED]The genetic structure of Aedes aegypti is driven by boat traffic

Discussion
Our results support the hypothesis that Ae. aegypti long-distance dispersal is driven by human activity, with particular emphasis on boat traffic. Low F ST values between Iquitos and surrounding towns in addition to our Structure results (i.e, the observation that no single population was exclusive to a single site) suggest mosquito mobility between human settlements. The pattern of mosquito gene flow also mirrors human transportation patterns: F ST values for a regional port (Iquitos B) were overall lower than F ST values collected from the interior of Iquitos (Iquitos C and D), implying more population mixing between Iquitos port mosquito populations and surrounding towns than between Iquitos interior populations and surrounding towns. Barrio Florida populations, in contrast, were the most genetically isolated compared to other towns. Iquitos is the major regional transportation hub, with frequent visits of barges and other types of boats from surrounding areas, whereas Barrio Florida is a small town (pop: 800), receiving most of its traffic from vehicles unlikely to be infested with Ae. aegypti (e.g., speed boats and small water taxis) [24]. Some limitations should be noted when interpreting our data. It is possible that patterns of human traffic, and by extension, Ae. aegypti gene flow, have changed since mosquito samples were collected in 2008. More recent samples and better transportation data would be required to assess this possibility. Our limited data on human transport patterns also meant that we were only able to evaluate the Propagule Pressure Index by comparisons between Iquitos and surrounding towns. To that end, in some locations our sample sizes were low, particularly in Aucayo. Further, some scholars have argued that the allozyme methods to estimate population subdivision (e.g., F ST ) are coarse measures of gene flow, and in some instances lead to erroneous conclusions [55]. Still, the results we present lay the groundwork for future population genetics approaches that may shed additional light on Ae. aegypti dispersal throughout the riverine landscape characteristic of this region. Fine-scale landscape genetics methods may be of particular interest to identify barriers to and pathways that facilitate dispersal [56].
Our previous studies in the Peruvian Amazon have shown that Ae. aegypti spreads farther along rivers than terrestrial routes [18], that Ae. aegypti spread is facilitated by certain vehicle types, large barges in particular (especially in comparison with terrestrial vehicles) [24], and that oviposition regularly occurs on boats while in transit [57]. When considered together with [EXSCINDED]The genetic structure of Aedes aegypti is driven by boat traffic past findings, this exploration of Ae. aegypti population genetics provides additional evidence to support the hypothesis that boats are major drivers of the passive transport of Ae. aegypti mosquitoes over long distances (between cities > 5 km apart). In other regions, fluvial traffic has also been implicated as a driver of Ae. aegypti spread [58], but the present study is the first to illustrate the importance of boat traffic in the Peruvian Amazon specifically.
Consistent with other findings from Peru [59], our comparisons of Ae. aegypti genetic variation and geographic distance did not reveal a pattern of isolation by distance. At this spatial scale (sampling locations < 100km apart), network distance might be the best predictor of genetic mixing. Although both field ecology [16,17] and population genetics studies [32,60,61] in the Americas have pointed to human transportation networks as a major driver of Ae. aegypti spread, to our knowledge, none of these studies linked ecological evidence with transportation and population genetics data (e.g., mosquito abundance collections, and data characterizing the number of vehicle trips between locations). The novel Propagule Pressure Index, PrPI, integrates these data by taking into account heterogeneity of Ae. aegypti infestation by vehicle type, in addition to the frequency of traffic between two locations. We were unable to statistically test for isolation by distance using PrPI because of limited sample sizes, however, the isolation by distance plot (Fig 2) is suggestive of a negative correlation between genetic and network distance. In other words, populations closer to one another in network space (as measured by PrPI) are likely to be more genetically related. Future studies with larger sample sizes are required to evaluate its utility.
Within Iquitos, Ae. aegypti genetic admixture of port mosquitoes and recently urbanized mosquitoes were similar, while interior mosquitoes had a different pattern of population membership. Within cities, mosquito population structure may be a function of availability of immature habitat, mosquito movement, or vector control programs that have differentially impacted various neighborhoods within the city. Perhaps, for example, insecticide application focused in the interior parts of the city act as barriers to population mixing, ultimately resulting in differing patterns of population group membership. This vector control explanation has been proposed in other settings [62].
Indeed, Ae. aegypti population structure in the Amazon has many other implications for dengue control programs. Our data suggest that Ae. aegypti are continually introduced to smaller surrounding communities from locations such as Iquitos that have stable and spatially diverse population structure. For example, the spread of insecticide resistance or vector competence genes would occur but could potentially lead to significant spatial heterogeneity. Significant spatial variation in vector competence has been reported previously [63], so understanding gene flow is critical for potential heterogeneity in transmission patterns. For new technologies that rely on population replacement (e.g., Wolbachia), it is possible that releases in large cities would be sufficient for spread to outlying areas.
An extension of our results is to ask, how does Ae. aegypti gene flow impact arboviral transmission dynamics? Dengue in the Peruvian Amazon has been characterized by transmission in large population centers with most cases observed from outlying communities assumed to be associated with frequent travel to these commercial hubs. The distribution of Ae. aegypti is driven by urbanization along highways, but infestation in communities along rivers is far less consistent spatially [18]. We hypothesize that establishment of Ae. aegypti in these communities is driven in part by repeated introductions as well as by local ecological characteristics such as water storage, management of containers, and microclimate. Over the last 20 years Ae. aegypti has spread to smaller communities, with significant implications for vector control programs which are already stretched to cover Iquitos and other larger cities. Dengue outbreaks are known to occur in several of the small riverine towns near Iquitos. Although these outbreaks are typically small and fleeting, the resources necessary for sustained control of Ae. aegypti (source reduction and larval control) are limited. Clearly, transmission in these smaller communities could lead to infections in travelers (and mosquitoes) that could then carry virus to other locations, thus serving as a mechanism of viral spread at broader scales.
Supporting information S1 Fig. DeltaK analysis of Structure results. The DeltaK statistic captures the degree of change between the log probability of the data for consecutive values of K populations. Results showed the most likely number of groups to be K = 3 genetic clusters. (TIFF) S1 Table. Pairwise distances between towns (Euclidean distance, fluvial path distance, shortest path distance, and PrPI). The number of visits between site pairs for large and medium barges were weighted according to the probability of their infestation (the columns Large Barges Weighted, Medium Barges Weighted). Our previous research demonstrated a 71% probability of infestation among large barges, 35% probability of infestation among medium barges, and 12.5% infestation among combis (taxis). PrPI is calculated by summing the number of visits by large and medium barges and combis, weighted for the probability of Ae. aegypti infestation. (DOCX) S2