Evidence for serial founder events during the colonization of North America by the yellow fever mosquito, Aedes aegypti

Abstract The Aedes aegypti mosquito first invaded the Americas about 500 years ago and today is a widely distributed invasive species and the primary vector for viruses causing dengue, chikungunya, Zika, and yellow fever. Here, we test the hypothesis that the North American colonization by Ae. aegypti occurred via a series of founder events. We present findings on genetic diversity, structure, and demographic history using data from 70 Ae. aegypti populations in North America that were genotyped at 12 microsatellite loci and/or ~20,000 single nucleotide polymorphisms, the largest genetic study of the region to date. We find evidence consistent with colonization driven by serial founder effect (SFE), with Florida as the putative source for a series of westward invasions. This scenario was supported by (1) a decrease in the genetic diversity of Ae. aegypti populations moving west, (2) a correlation between pairwise genetic and geographic distances, and (3) demographic analysis based on allele frequencies. A few Ae. aegypti populations on the west coast do not follow the general trend, likely due to a recent and distinct invasion history. We argue that SFE provides a helpful albeit simplified model for the movement of Ae. aegypti across North America, with outlier populations warranting further investigation.

, with a rate of invasion that is expected to accelerate with climate change (Iwamura et al., 2020). The global distribution of this vector enabled the recent outbreaks of Zika and chikungunya in the Americas (Carlson et al., 2016;Leparc-Goffart et al., 2014), reemergence of yellow fever in Africa and South America (Hamlet et al., 2018), and a dramatic spread and increase of dengue cases around the world (Brady & Hay, 2020). Its widespread presence in Mexico and the Caribbean perpetuates endemic dengue, and its presence in the southern United States poses a public health threat.
Aedes aegypti likely arrived in the Americas in the 17th century aboard slave ships from Africa, where it rapidly spread throughout the continent, as evidenced by outbreaks of yellow fever and dengue; see the comprehensive review . Disease outbreaks ranged from the Caribbean to the North Atlantic and the southeast United States by the early 1800s (Carrigan, 1959;Moreno-Madriñán & Turell, 2018). Dengue made its way into the central-south United States in the 1850s (Chandler, 1956), and by the 1930s Ae. aegypti spanned much of Texas and coastal Mexico (Slosek, 1986). Ae. aegypti became more prevalent and widespread in New Mexico around 1994 (Merrill et al., 2005) and in California starting around 2013 (Metzger et al., 2017). Taken together, these lines of evidence paint a picture of the westward migration of the mosquito species.
When a small number of founders emigrate from a larger population, the new population will generally display a reduction in genetic variation relative to the original population that can last for many generations-a phenomenon called "founder effect" (Nei et al., 1975). The serial founder effect (SFE) (spread occurring through a series of these founder events) has been famously invoked to describe the movement of Homo sapiens out of Africa (Henn et al., 2012;Ramachandran et al., 2005). It has also been used to explain the global spread of the malaria parasite that accompanied humans out of Africa and the expansion of monarch butterflies across the Pacific (Pierce et al., 2014). Although numerous studies have examined Ae. aegypti invasions, migration, and structure using populations genetics from around the globe-including North America (Gloria-Soria et al., 2016;Kotsakiozi et al., 2018;Pless et al., 2017)none have explicitly tested the hypothesis of colonization by SFE.
Here, we asked whether SFE explains the westward spread of Ae. aegypti across southern North America. These migrations were likely accomplished by a combination of active dispersal and passive human-mediated transport (e.g., trucks and ships) (Fonzi et al., 2015; Goncalves da Silva et al., 2012;Guagliardo et al., 2014;Medley et al., 2015). Given an active average lifetime dispersal of <200 m for Ae.
aegypti (Honorio et al., 2003;Jasper et al., 2019;Reiter, 2007;Russell et al., 2005), we predict that short-distance invasions are more common for both active and passive dispersal due to a higher number of introduction events and higher propagule pressure (Sakai et al., 2001), which is likely to result in colonization via SFE.
To test if Ae. aegypti expansion across North America is consistent with a SFE moving westward, we evaluated our data against four expectations: 1. Source populations will have the highest number of alleles, with genetic diversity decreasing in proportion to the distance from the source (Hunley et al., 2012).
2. There will be a positive relationship between distance among groups and their genetic differentiation (Ramachandran et al., 2005).
3. Following each founding event, the daughter group will carry a subset of variation from the parental group (Ramachandran et al., 2005).
4. Demographic inference modeling will support SFE over other possible scenarios.
We inferred genetic structure across the region and tested these predictions using data from 12 microsatellite loci and single nucleotide polymorphism (SNP) array data from 70 North American Ae. aegypti populations. Establishing the extent to which Ae. aegypti has spread via SFE is important for preventing and detecting future invasions, as well as modifying vector control in response to the movement of pesticide-resistant alleles.

| Mosquito collection
Our analysis includes 70 Ae. aegypti populations across continental North America and the eastern Caribbean ( Figure 1, Table 1 and   Table S1). Microsatellite and SNP genotypes for most populations in this study have been reported in Evans et al. (2015), 2007). The SNP dataset includes 373 individuals from 36 populations genotyped at ~20,000 SNPs with the Axiom_aegypti array (Evans et al., 2015). Microsatellites are appropriate for unbiased genetic diversity estimates and demographic inference analysis due to their multiallelic, highly polymorphic nature, and complete allele frequency spectra. Additionally, their low cost facilitates the generation of large global and local reference databases. In contrast, the larger number of markers from the Axiom_aegypti array may provide higher fine-scale resolution for population structure, particularly when populations are of recent origin or when gene flow is significant . However, the ascertainment bias in the design of the SNP array is likely to affect measurements of genetic diversity and may also impact demographic analysis, which uses allele diversity spectrum to evaluate the likelihood of different scenarios. Therefore, the two sets of genetic markers provide complementary advantages and are useful for different analyses.
Whole genomic DNA was extracted from individual mosquitoes using the Qiagen DNeasy Blood and Tissue kit according to manufacturer instructions, including the optional RNAse A step. Individuals were genotyped for 12 microsatellites as in Gloria-Soria et al. (2016).

F I G U R E 1 (a) Map showing locations for each
Aedes aegypti sampling site included in this study, colored by regional group. The inset on the bottom left shows southern California. Mean allelic richness estimated from microsatellites by rarefaction of all sites within each region is shown in the boxes. (b) STRUCTURE plot of the complete microsatellite dataset with K = 4 number of clusters. Each column is an individual, and the heights of the color bars represent the proportion of ancestry that came from each of the four inferred clusters (yellow, light blue, red, and dark blue). The outlier (red) in southern California is Exeter County. The longitude of each site is plotted against its mean observed heterozygosity (c), expected heterozygosity (d), and allelic richness estimated by rarefaction using microsatellite data (e). Colors in the map (a), above the STRUCTURE plot (b), and in the diversity plots (c-e) are consistent and correspond to regional groups Four loci (A1, B2, B3, and A9) are trinucleotide repeats, and eight (AC2, CT2, AG2, AC4, AC1, AC5, AG1, and AG4) are di-nucleotide repeats. Any individuals that were genotyped at fewer than 10 loci were excluded from the analysis. We include only populations with a minimum of 10 individuals and arbitrarily selected 40 individuals from sites with more than 40 samples to control for unequal or low sample size (Puechmaille, 2016).
Individuals were genotyped using Axiom_aegypti, a highthroughput genotyping chip that has 50,000 probes (Evans et al., 2015). Genotyping was conducted by the Functional Genomics Core at the University of North Carolina, Chapel Hill. To prune the SNP dataset, we first excluded 2166 SNPs that failed a test of Mendelian inheritance (Evans et al., 2015). Since some analyses can be confounded by SNPs in linkage disequilibrium (Alexander et al., 2009), we excluded tightly linked SNPs with the plink command "--indep 75 kb 50 2" (Purcell, 2016;Purcell et al., 2007). We also excluded any SNPs that were genotyped in <90% of the individuals and those with a minor allele frequency of <5%, resulting in 20,003 SNPs remaining for analysis.
All microsatellite data are available in Dataset S1, and all SNP data are available in Dataset S2. Additionally, the data can be accessed on Dryad (https://doi.org/10.5061/dryad.5x69p 8d5j) and VectorBase (VBP0000801).

| Inferring geographic regions based on genetic structure
All microsatellite loci were tested for within-population deviations from Hardy-Weinberg equilibrium and for linkage disequilibrium among loci pairs using the R package Genepop v. 1.1.4. with 5000 dememorizations, 500 batches, and 5000 iterations per batch for both tests (Raymond & Rousset, 1995). To correct for multiple testing, a Bonferroni correction was applied at the .05 α level of significance.
To establish regional groupings for subsequent analysis, we examined population structure using a number of methods. We performed principal component analysis (PCA) using the R pack- Note: Population number corresponding to Figure 1, population name, population abbreviation, region name, year sampled, the sample size for microsatellite data (N 1 ), and sample size for SNP data (N 2 ). New microsatellite data are indicated with an asterisk (*) and new SNP data are indicated with a cross ( + ).
Abbreviations: DR, Dominican Republic; MEX, Mexico; USA, United States of America; USVI, United States Virgin Islands.

TA B L E 1 (Continued)
We conducted 20 independent runs of STRUCTURE v. 2.3.4 for K = 1-10 (Pritchard et al., 2000) using the microsatellite data and 600,000 generations, with the first 100,000 discarded as burn-in.
We visualized the STRUCTURE results using the programs Clumpak and DISTRUCT v.1.1 (Kopelman et al., 2015;Rosenberg, 2004), and we inferred the optimal value of K using relevant guidelines (Cullingham et al., 2020;Earl, 2012;Evanno et al., 2005). We repeated these analyses for the Caribbean microsatellite dataset (40 independent runs of STRUCTURE for K = 1-6), and for the Central and Southwest microsatellite datasets combined (20 independent runs of STRUCTURE for K = 1-10) to further explore the regional genetic structure of the two regions that include new populations.
Additionally, we evaluated genetic structure within the Caribbean using a multivariate approach, Discriminant Analysis of Principal Components (DAPC), on the microsatellite data using the Adegenet package (Jombart, 2008).
Based on the results from our genetic clustering analysis, we grouped the samples into six regions for further analyses (

| Prediction 1: Genetic diversity decreases toward the west
Observed heterozygosity (H O ), expected heterozygosity (H E ), the inbreeding coefficient (F IS ), and a number of private alleles (alleles found in no other population) were calculated from the microsatellite dataset for each population using GenAlEx v. 6.51 (Peakall & Smouse, 2006), and allelic richness was estimated by rarefaction (N = 30) using the software HP-Rare v. 1.0 (Kalinowski, 2005). These measurements were not calculated for the SNP dataset due to ascertainment bias in the design of the SNP array (Evans et al., 2015). We assessed regional genetic diversity by averaging the genetic diversity metrics across all the sites within each region. Since regional allelic richness and the number of private alleles could be biased by the number of individuals in a region, we calculated the genetic diversity measures a second time after combing all individuals within each region and then randomly resampling them so that each region had the same number of individuals.

| Prediction 2: Positive relationship between geographic and genetic distance
Using both the microsatellite and SNP datasets, we calculated pairwise F ST and evaluated significance with 1000 permutations using Arlequin 3.5 (Excoffier et al., 2005). We then tested for a relationship between pairwise genetic distance Because clustering and assignment methods may mistake continuous processes (e.g., isolation by distance, in which there is a positive relationship between geographic distance and gene flow between populations) for discrete processes, we implemented a method called conStruct v. 1.0.4 (Bradburd et al., 2018) using the SNP array data. This model-based clustering method uses isolation by distance when possible to explain genetic variation. We ran three independent runs of the program for K = 1-4 with 1000 iterations and the spatial model setting and another three independent runs with the non-spatial model setting. We assessed posterior probability, Markov chain Monte Carlo (MCMC) performance, and layer contributions to compare the independent runs and identify a putative optimal number of layers for each run. The program performs best when there are more loci than the number of samples, so it was not suitable to run with the microsatellite data.

| Prediction 3: Daughter groups nested within parental group
To test the hypothesis that daughter groups contain a subset of the allelic diversity from parent groups, we created the presence/ absence matrices representing all alleles for each of the 12 microsatellites, where regional sites were columns, and the different alleles were rows. Microsatellites were used rather than SNPs because the latter are biallelic, and thus do not have the broad allele frequency spectra required for this type of analysis. To maximize the chance of finding a pattern if one existed, we focused on three regions that were most likely to display this pattern: Southeast, Central, and Southwest. To control for uneven sample size across different populations, we combined individuals within each region and resampled, so each region was represented by the same number of individuals before analysis. We calculated the nestedness metric based on the overlap and decreasing fill (NODF) for each matrix, in which higher scores indicate greater amounts of nestedness (Morrison, 2013), using the R package "RInSp" (Zaccarelli et al., 2013). To compare our matrices and their NODF scores with the null expectation, we created five control model matrices by shuffling the elements in each row using the R package "picante," specifically the randomizeMatrix tool, with 1000 iterations and the null model set as "richness" (Kembel et al., 2010).

| Prediction 4: Demographic inference modeling supports SFE
We performed demographic history analysis using DIYABC-RF Priors were set to be as wide as possible within reason given known colonization events (Table S2). Each resulting dataset is then summarized with a set of descriptive statistics. These summary statis- Specifically, we compared SFE scenario-in which Southern CA split from Southwest, which split from Central, which split from the Southeast-against a scenario where each region diverged from the Southeast independently ( Figure 2). Each region was represented by 158 randomly selected samples, and we performed the analysis in duplicate, using a second random draw from each region. We excluded Exeter from Southern California since it is a known outlier that likely has a different origin than the other sites in Southern California (see Results). To ensure that pooling across populations did not bias the outcome, we also repeated the analysis twice replacing the regions with single, arbitrarily selected populations from each region (Table S2). We ran DIYABC-RF using a training set of 200,000 simulated datasets. The constructed RF for both model choice and parameter estimation used 1000 trees, and the best scenario was selected based on linear discriminant analysis and partial least squares regression analysis, in each case. Full parameters and priors are listed in Table S2.

| Inferring geographic regions based on genetic structure
We first confirmed that our microsatellites behave as independent, To test our hypothesis that Ae. aegypti primarily colonized North America through a series of westward founder effects, it is useful to represent our 70 populations as major regions defined by their geography and population structure. We explored genetic population structure using both microsatellites and SNP markers and various methods. Using Bayesian clustering with the CV error F I G U R E 2 Scenarios tested by DIYABC-RF using microsatellite data. The two trees depict scenarios of Ae. aegypti colonization of continental North America. Scenario 1 shows independent invasions from the Southeast with bottlenecks after each lineage division and Scenario 2 shows a pattern of serial founder effect moving west from the Southeast, also with bottlenecks. Posterior probabilities (PP) are shown for each scenario. Time 0 corresponds to the present (when the samples were collected). Colors for each lineage correspond to Figure  1, and changes in effective population size (bottlenecks) are shown in gray. Posterior probabilities and divergence times correspond to a single run of the model using random draws from each region (Table S2b); results are similar for the other three independent runs (Table S2) method (Evanno et al., 2005) for the microsatellite data, we found two primary clusters which generally split the west half of North

America from the east half, with a large transition zone through New
Mexico and Texas ( Figure S1a). At K = 4, the four ancestries correspond generally to 1: Northern California + Central, 2: Southwest, 3: Southeast, and 4: Caribbean (Figure 1b). Southern CA has sites that are outliers (e.g., Exeter) or appear to show admixture (e.g., Mission Viejo), consistent with the complex, recent invasions known to characterize the region (Pless et al., 2020).

| Prediction 1: Genetic diversity decreases toward the west
The microsatellite data showed that, overall, the highest genetic diversity was found in the Southeast, followed by the Caribbean, Central, Southwest, Northern CA, and Southern CA (R 2 for longitude vs. allelic richness = 0.32, p < 10 −8 ) (Figure 1, Table 2). After combining and standardizing the number of individuals per region, we found that Southeast and Southern California had the highest number of private alleles (alleles found in no other region), followed by the Caribbean and Central (Table 2). This is consistent with the Southeast having the oldest population and Southern California being founded multiple times from different regions. There was a stronger correlation between heterozygosity versus distance from Florida (R 2 = .53, p < 10 −11 ) than between heterozygosity and distance from the Caribbean (R 2 = .32, p < 10 −11 ).

| Prediction 2: Positive relationship between geographic and genetic distance
Considering the microsatellite data, the F ST (a measure of genetic differentiation due to genetic structure) values between all population pairs (N = 1953) were 0.13 ± 0.070 (mean ± SD), and all values were significantly different than zero at the 0.05 level after correction for multiple tests (Table S4). Geographic distance was a significant predictor of genetic distance (calculated as To investigate the effect of continuous processes, namely isolation by distance, on clustering, we ran conStruct for K = 1-4 using the SNP dataset. We focus on K = 2 because this model consistently produced well-behaved MCMC posterior probability plots, and the first two layers contributed most of the variance in the majority of the runs (Table S6). The three independent non-spatial models produced very consistent results at K = 2, showing a clear structure, especially between Florida + Exeter and California ( Figure 5 and Figure S4a). After controlling for the effect of isolation by distance, the three spatial runs showed significantly less genetic structure ( Figure 5 and Figure S4b), although one run continued to show Exeter as a clear outlier ( Figure S4c). In all cases, the posterior probability across MCMC runs was higher for the spatial models than the non-spatial models. The non-spatial model showed additional structure at K = 3 and K = 4 (e.g., differences between the Caribbean and continental North America), but again the spatial model at these higher K levels showed less structure.

| Prediction 3: Daughter groups nested within parental group
To evaluate if each putative daughter group contained a subset of the diversity in the presumed parent group, we calculated NODF for each microsatellite allele for the three regions where we expected the larger effect: Southeast, Central, and the Southwest. However, nine loci had <4 alleles, an insufficient number to detect nestedness patterns (Table 3). Of the three remaining loci, one (AG2) showed evidence of nestedness (NODF score >1SD above the neutral NODF scores), and the other two were within 1SD of the neutral model (Table 3). Additional analyses using more regions or individual populations rather than regions did not show significant nestedness.

| Prediction 4: Demographic inference modeling supports SFE
Using approximate Bayesian computation combined with supervised machine learning, we compared two invasion scenarios to test whether the data were more consistent with SFE or independent invasions from the Southeast (Figure 2). SFE had a higher posterior probability (PP > 0.922 across all four independent runs) than independent invasions from the Southeast (PP < 0.07 across all four independent runs). Under the SFE model and assuming ten genera- tions/year, we estimated that Southern CA split from Southwest ~39-66 years ago, Southwest split from Central ~166-207 years ago, and Central split from Southeast ~322-345 years ago (Table S2).
We ran the model twice using randomly selected samples from each region and twice using single populations from each region; results were similar across the four runs (Table S2).

| Regional genetic structure
In this study, we tested whether SFE described the westward spread of Ae. aegypti across North America by evaluating the genetic diversity and population structure from a rich dataset of 70 populations.
We first established relevant regions for analysis based on genetic structure; these are the hypothesized "stepping-stone" regions for the westward spread of Ae. aegypti. We identified six genetic regions F I G U R E 3 Population structure in North America generated by fastStructure using the SNP dataset. Populations are arranged according to their longitude (west on the left). Each column is an individual, each color represents an inferred group, and the height of the color bar shows the proportion of ancestry that came from the K inferred groups (K = 2 on the top panel and K = 8 on the bottom panel) based on both microsatellites and SNPs that are consistent with previous work (Gloria-Soria et al., 2016;Kotsakiozi et al., 2018). Broadly, we detected two primary clusters in the microsatellite dataset which split the eastern populations from the western ones with a large admixture zone through the Southwest and Central regions (Figure 2a).

TA B L E 2 Genetic diversity of each region calculated with microsatellites
At higher levels of K, additional structure emerged (Figure 2), and we found genetic differentiation by geography with some outliers, specifically Exeter CA, which clusters with Florida and likely represents the third invasion into the state ( Figure S1). Southern California has unusual genetic patterns (e.g., high genetic differentiation and a high number of private alleles despite low diversity), as discussed elsewhere (Pless et al., 2020). The stability of the regions we identi- the islands as previously reported (Wallis et al., 1984). Although Ae.
aegypti was likely established in the Caribbean and Florida at similar times in history (1600-1700s), the lower genetic diversity and higher amount of genetic structure in the Caribbean may reflect the accumulation of genetic drift due to smaller population sizes and lower migration between the islands than between sites in Florida and/or population bottlenecks from eradication efforts on the islands (Dick et al., 2012;Sherpa et al., 2018).

| Genetic predictions of SFE
The genetic data supported three of our predictions of SFE (#1, #2, and #4), with the Southeast as the putative source for westward spread across North America. The highest numbers of private alleles were found in Florida and the Caribbean, and genetic diversity largely decreased in populations westward (Tables 1 and 2, Figure 1).
Distance to Florida was strongly correlated with heterozygosity; distance to the Caribbean was correlated to heterozygosity, but the signal was weaker. Florida likely harbors the oldest and most stable populations of Ae. aegypti in continental North America today because (1) southern Florida is hospitable year-round for Ae. aegypti (Johnson et al., 2017), (2) vector control never fully eliminated Ae.
In line with our second expectation of SFE, geographic distance was a strong predictor of genetic distance, especially after removing recent invasions on the west coast (including the highly differentiated southern California sites, and the northern California sites, which were likely founded by long-distance human-mediated movement)   (Figure 4). We also found that some of the genetic structures in North America could be explained by a continuous process like isolation by distance (Bradburd et al., 2018).
Accounting for genetic variance due to isolation by distance resulted in less genetic clustering ( Figure 5 and Figure S4).
Our third expectation, that daughter populations would contain a subset of the diversity contained by their parental group, could not be properly tested. The overall low genetic diversity of Ae. aegypti in the Americas, relative to that of its native Africa (Gloria-Soria et al.,

2016), resulted in a limited number of alleles per microsatellite and
frequent presence of rare alleles that prevented the detection of diversity nestedness in the dataset. Exploring different regional categories and using higher sample sizes could help clarify these results.
We do find support for the SFE model in the demographic analysis (prediction 4). Specifically, westward spread via a series of founder effects was better supported than independent invasions from the Southeast.
To the best of our knowledge, the SFE model has not been tested for any mosquito species, and has only been tested on a few insects more generally, see Pierce et al. (2014). We argue this model is useful in thinking about the expansion of Ae. aegypti out of Africa and across continents such as North America. Like isolation by distance, SFE can serve as a null model for Ae. aegypti genetic structure in North America, and outliers from the model warrant additional research and explanation (e.g., California being founded by multiple invasions, some from across the country). Our study also has some limitations that should be addressed in future work. While the linear decline of heterozygosity is indicative of SFE, it can also be explained by extensive admixture from an early branching lineage into later branching lineages (Pickrell & Reich, 2014). The increase in genetic distance with an increase in geographic distance can be caused by the budding pattern of migration under SFE or solely by isolation by distance (Wright, 1943). More extensive sequencing and modeling would be useful for testing these different demographic scenarios.
Clearly, Ae. aegypti can and does make stochastic long-distance jumps, for example, from the U.S. southeast to the Netherlands (Brown et al., 2011) and to California (Gloria-Soria et al., 2014).
However, the demonstrated relevance of SFE and of environmental Note: Asterisk (*) indicates nestedness metric (NODF) significantly greater than neutral score.
Abbreviations: Microsatellite, each microsatellite locus; NODF, nestedness metric based on the overlap and decreasing fill; Neutral ± SD, generated by randomly shuffling the rows of the nestedness matrix; No. alleles total, number of alleles for given microsatellite; No. alleles not fixed, number of alleles that vary across the three regions.

ACK N OWLED G M ENTS
We thank our sample contributors: Kathryn Hanley and the NMSU- Services.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genetic data (microsatellite calls and SNP array dataset) are included in the supplementary materials and are archived at Dryad (https://doi.org/10.5061/dryad.5x69p 8d5j) and VectorBase (VBP0000801).