Evidence of Local Extinction and Reintroduction of Aedes aegypti in Exeter, California

Established populations of Aedes aegypti, a mosquito vector of multiple major arthropod-borne viruses, were first found in three California (CA) cities in 2013. From 2013 to April 2021, Ae. aegypti thwarted almost all control efforts to stop its spread and expanded its range to 308 cities, including Exeter, in 22 counties in CA. Population genomic analyses have suggested that multiple genetically distinct Ae. aegypti populations were introduced into CA. However Ae. aegypti collected for the first time in 2014 in Exeter, appeared to be different from three major genetic clusters found elsewhere in CA. Due to intense control efforts by the Delta Vector Control District (DVCD), Ae. aegypti was thought to have been eliminated from Exeter in 2015. Unfortunately, it was recollected in 2018. It was not clear if the reemergence of Ae. aegypti in Exeter was derived from the bottlenecked remnants of the original 2014 Exeter population or from an independent invasion from a different population derived from surrounding areas. The goal of this work was to determine which of these scenarios occurred (recovery after bottleneck or reintroduction after elimination) and if elimination and reintroduction occurred to identify the origin of the invading population using a population genomic approach. Our results support the reintroduction after elimination hypothesis. The source of reintroduction, however, was unexpectedly from the southern CA cluster rather than from other two geographically closer central CA genetic clusters. We also conducted a knockdown resistance mutation profile, which showed Exeter 2014 had the lowest level of resistant alleles compared to the other populations, could have contributed towards DVCD’s ability to locally eliminate Ae. aegypti in 2014.


INTRODUCTION
Aedes aegypti serves as a major vector of four human diseasecausing viruses, including yellow fever, dengue, chikungunya, and Zika viruses, posing a major threat to public health. Records indicate this species became established in the southeastern United States of America between the 15-18th centuries (1) and then spread throughout the east coast and southern states (2). California (CA) had remained free from Ae. aegypti until the summer of 2013 (3,4) when the species were first detected in Menlo Park, Clovis, and Madera (4). Since then, this species has expanded its range to 308 cities in 22 counties, as of April, 2021 (5).
Attempts to locally eliminate and even control this highly invasive species has proven to be extremely challenging. For example, the Consolidated Mosquito Abatement District (CMAD) implemented labor intensive integrated vector control management in the city of Clovis where Ae. aegypti were first detected in 2013. Their efforts involved extensive public education, thorough property inspections, sanitation, insecticide treatment at larval sources, and residual barrier spraying with pyrethroids. Despite these efforts, Ae. aegypti successfully overwintered and continued to persist in Clovis (6).
Delta Vector Control District (DVCD), which covers northern Tulare County, just south of the area covered by CMAD, first detected Ae. aegypti in Exeter in August of 2014. Following detection, the DVCD initiated an intensive control campaign involving thorough, routine property inspections and barrier applications, public education, breeding site treatment or removal, as well as hand and truck-mounted adulticide fogging through October. DVCD encountered significant pushback from the public due to the regular adulticide applications, mandatory property inspections, and confiscation of container habitats. However, they reported that residents followed instructions and were able to control mosquito breeding. Exeter remained free of Ae. aegypti from the beginning of 2015 through the summer of 2017. DVCD detected Ae. aegypti again in 2017 in the neighboring cities of Visalia and Farmersville and, in 2018, at multiple sites in Exeter. In response to this detection, the district attempted to mount a response similar to 2014. Despite frequent, mandatory property inspections, breeding site elimination and sanitation, and barrier spraying, the infestation persisted.
Population genetic/genomic analyses suggest that multiple genetically distinct Ae. aegypti populations were introduced in CA (7, 8). The Ae. aegypti populations in California could be largely grouped into three major genetic clusters (8). One cluster includes samples from Fresno, Madera, and Menlo Park in Central CA (8). Another cluster includes samples from Clovis in central CA adjacent to Fresno. The third cluster includes all southern CA samples, as well as 2014 Exeter samples (8). Further clustering revealed that 2014 Exeter samples are similar to an Ae. aegypti population from Florida rather than the populations from southern CA. However, it is important to note that the data from Lee et al. (8) does not provide definitive evidence that Exeter Ae. aegypti were introduced from Florida because the cluster analysis included only Florida samples for comparison. Aedes aegypti from southeastern states such as Louisiana and Texas are genetically similar to Florida populations based on microsatellite analysis (3).
The genetic dissimilarity of 2014 Exeter to other California Ae. aegypti provided an opportunity to investigate the hypothesis that (1) the Exeter mosquito population went through a severe bottleneck due to intensive insecticide spraying and control followed by a recovery after a couple of years (bottleneck and recovery) or (2) DVCD successfully eliminated its initial introduction but later new mosquitoes migrated or were reintroduced to Exeter from other cities (local extinction and reintroduction). If bottleneck and recovery occurred, then we expected 2018 Exeter samples to have similar genetic profile as 2014 Exeter samples. If localized extinction and reintroduction occurred, then we expected 2018 Exeter samples to have a similar genetic profile to any of the other three CA genetic clusters. To investigate this hypothesis, we report the population genomic analysis of 243 Ae. aegypti from California, Arizona, Florida, and Mexico in this paper.

Sample Collection
Genome data was obtained from deposited sequences available for specimens originating from earlier collections in Clovis, Fresno, Madera, Menlo Park, 2014 Exeter, East Los Angeles, San Diego (CA), Vero Beach, and Key West, Florida (FL) -NCBI BioProject PRJNA385349 (8,9). Further samples from St. Augustine (FL), Naples (FL), St. Lucie (FL), Arizona, and Mexico were collected as adults using BG Sentinel traps and the remaining samples as eggs in ovicups, which were then reared to the adult stage in the laboratory and stored in >70% alcohol prior to DNA extraction. The global positioning system (GPS) coordinates from where each site specimens originated from as well as the number of samples sequenced and genotyped are provided in Table 1. The sample locations map is provided in Figure 1.

Genome Sequencing
Protocols for sequenced specimens available in the NCBI Bioproject database are described in Lee et al. (8) and Schmidt et al. (9). For the other specimens, DNA was extracted from individual mosquitoes using a magnetic-bead based DNA extraction protocol described in Chen et al. (10)

SNP Genotyping
Eggs were collected from the field in Clovis and Greater LA and reared in the lab under existing protocols (11). Adult collections from Exeter were obtained from Delta Vector Control District. Sixty individuals from Clovis, 67 individuals from Greater LA, 24 individuals from Exeter 2014 and 24 individuals from Exeter 2018 were collected ( Table 1, N KDR ). DNA was extracted using the Zymo Quick-DNA/RNA Mini Prep Kit (#D7001) using the protocol for Solid Tissue. DNA quality and quantity were determined using a Qubit instrument (Life Technologies) and approximately 4 ng/µL was extracted from each individual.  DNA was submitted to the UC Davis Veterinary Genetics Laboratory for SNP genotyping using the iPLEX MassARRAY analysis following the protocol described in Lee et al. (12). Thirty-seven SNPs were selected from the published whole genome sequences (8), 29 of which were chosen due to their association with specific genetic clusters and eight within the Voltage Gated Sodium Channel gene. Data for 5 knockdown resistance (kdr) SNPs (F1534C, V410L, S723T, I915K, and V1016I) were included from (13). A full list of SNPs and primers can be found in Table S1. Populations were clustered using a Principal Component Analysis (PCA) performed in R v4.0.5, removing individuals with no calls for any SNPs.

Genome Sequence Data Analysis
Raw reads were trimmed using fastp (14) version 0.20.1. Trimmed reads were mapped to the Ae13CLOV028MT (Genbank ID: MH348176) first using BWA-MEM (15) version 0.7.15 following recommendation from Schmidt et al. (16) to minimize the impact of mitochondrial reads mapping to the nuclear genome due to presence of pseudogenes (17). Unmapped and mate-is-unmapped reads from mitogenome mapping were filtered using sambamba, converted to fastq files using samtools version 1.12, and mapped to AaegL5 reference genome (18) using BWA-MEM (15)  The repeat regions were soft-masked in the AaegL5 reference genome and SNPs in these regions were excluded from analysis. Only biallelic SNPs with a minimum of 6X coverager were used for further analysis. A missing data threshold of 10% was used to filter SNPs. Hudson FST (21) and Principal Component Analysis (PCA) analyses was done in Python version 3.6.6 using the scikitallel module version 1.2.0 (22). A phylogenetic tree based on the SNPs was constructed using the neighbor-joining algorithm as implemented in PHYLIP (23) version 3.696. Bayesian clustering method implemented in ADMIXTURE v1.3.0 (24) was used to estimate ancestry components for each individual. For this analysis, a total of 10 iterations were performed for values of K clusters from 1 to 10 with no prior population assignment. The results for each K were compiled using the online version of CLUMPAK and plotted in R. Windowed population genetic statistics such as F ST and nucleotide diversity (p) were calculated using scikit-allel library with a 1Mbp window size and half-step overlapping windows.

RESULTS
Thirty-eight Ae. aegypti genomes were sequenced, originating from Florida (N=24), Arizona (N=4), Exeter 2018 (N=4), and Clovis (N=6). These genomes were analyzed together with 30 sequenced genomes used in other publications (N=68) ( Table 1, N WGS ) (8,9). Each sample was sequenced with a mean nuclear genome coverage of 11.9X per sample (Table S2; Figure 2). One Exeter 2018 mosquito appears to have a signature of an admixture with the Clovis cluster. This particular sample corresponds to the outlier individual in the PCA that occupied a space between the southern CA and central CA clusters (Figure 1). The likelihood values increase continually with higher K ( Figure S1) indicating that further substructure may be present beyond the K=5 clusters presented in Figure 2.
Exeter 2014 and 2018 were the most distant population pairs (F ST = 0.152; Figure 3 and Table S3). Consistent with Admixture results, East Los Angeles, Arizona, Mexico, and Florida form a closely related group with F ST < 0.05 (Table S3). Fresno, Madera, and Menlo Park also showed no differentiation with F ST = 0 (Figures 2, 3 and Table S3). The closest population pair with Exeter 2014 was Florida (F ST = 0.055) while the closest population pair with Exeter 2018 was East Los Angeles (F ST = 0.038).
Based on published genome sequences (8), we selected a set of biallelic SNPs most informative for separating individuals into each of the major three genetic clusters found in CA. We designed a multiplex SNP genotyping assay using the iPLEX MassARRAY system to screen 37 SNPs simultaneously for their genetic background (Table S1). Thirty-one of these SNPs differentiate Ae. aegypti into each of the four genetic clusters present in California and 6 of these SNPs are associated with permethrin resistance. We genotyped 2018 Clovis, Greater Los Angeles, Exeter 2014, and Exeter 2018 samples using this new multiplex SNP assay. Our PCA analysis of SNP genotypes (Figure 4) shows Clovis and Greater Los Angeles forming separate genetic clusters, although the boundary is not as clear as genome-wide SNP data shown in Figure 1. Consistent with genome-wide SNP data, Exeter 2018 has large overlap with southern CA samples with two samples potentially having mixed ancestry from the southern CA and Clovis clusters. Exeter 2014 Ae. aegypti clearly clustered separately from the Clovis 2018, Exeter 2018, and Greater Los Angeles Ae. aegypti mosquitoes.
We calculated the windowed F ST , nucleotide diversity (p), and change in p between 2018 and 2014 Exeter samples ( Figure 5) to identify any hotspots of divergence. Any F ST greater than 0.1 (red line in Figure 5A) indicates high genetic distance in the corresponding genomic region between two years. The genomic divergence between the two years appears to be genome-wide except in limited locations such as in the immediate vicinity of the centromeres. Overall elevated p was observed near the centromeres of Chromosome 1 and 3 ( Figure 5B). While 2018 Exeter samples had higher p in Chromosome 1, 2014 Exeter samples had higher p in Chromosome 3 ( Figure 5C).
We also analyzed 6 SNPs in the voltage-gated sodium channel gene, typically known as knockdown resistance gene (kdr) in mosquito literature (AAEL023266) located on Chromosome 3 using the method described in Mack et al. (13). Five of these 6 SNPs were previously used in Mack et al. (13). The SNP Q1853R (3:315931672) is a new addition to this study as it was only detected in the mosquitoes derived from the Exeter 2014 population. None of the other Californian populations possess this SNP and it appears to have disappeared with the eradication of that population. More than 60% of Exeter 2014 and 2018 Ae. aegypti had a V1016 resistant allele, while >90% of Clovis mosquitoes had this resistant allele present. Other nonsynonymous mutations in the kdr gene also occurred at the highest frequencies in Clovis. Except for allele F1534C, which was almost fixed in Clovis, alleles V410L, S723T and I915K were at much lower frequencies in Exeter 2014, 2018 and Greater LA Ae. aegypti ( Table 2).   Had the Exeter population experienced a severe reduction due to DVCD control efforts, followed by 3-4 years of recovery before detection again, then the 2018 Exeter individuals would be expected to be genetically similar to those from Exeter 2014. However, the 2018 Exeter samples form a genetic cluster separate from that of the 2014 Exeter cluster. Genetic differentiation between 2014 and 2018 Exeter samples appears genome-wide with the majority of genomic regions exceeding F ST > 0.1 ( Figure 5). This pattern is in contrast to the level of genetic differentiation we observed between Clovis 2013 and 2016 (8). Such a drastic change in genome structure would be extremely unlikely to develop within a few years simply by genetic drift or natural selection, adding additional evidence toward the local elimination and reintroduction hypothesis for Ae. aegypti population in Exeter, CA. The Exeter 2018 cluster more closely aligns with the samples from the Southwestern USA and Mexico indicating that 2018 samples likely resulted from a reintroduction to Exeter. The most likely source of reinvading individuals would intuitively come from neighboring cities like Clovis and Fresno, which experience high abundances of Ae. aegypti in the middle to late summer. However, the genetic similarity of the Exeter 2018 to southern CA Ae. aegypti suggests that the founding Exeter 2018 mosquito/ es were likely escapees traveling in a vehicle traveling from southern CA, demonstrating the propensity of Ae. aegypti for human-mediated dispersal across long distances. Our study is the first report of Ae. aegypti from a Central CA location sharing the majority of its genetic profile with Ae. aegypti from southern CA. Importantly, when Ae. aegypti were detected in 2018 in Exeter, they were also detected in several Visalia suburbs (unpublished data), and it remains to be learned whether the Ae. aegypti that are now present in Visalia and other towns neighboring Exeter are genetically similar or different to those from Exeter and whether they are migrants from the second Exeter invasion or from other neighboring central California locations. The success in eliminating the population detected in 2014 is likely due to the limited geographical scope of the population and the intensive control efforts that were implemented. However, the possibility of eliminating the new populations of Ae. aegypti in Exeter and suburbs of Visalia is unlikely as resources are stretched too thin to mount such an intense control effort across such an extensive area.
Southern CA Ae. aegypti had different frequencies of nonsynonymous SNPs or mutations in the VGSC gene than in Central CA populations, including the ubiquitous V410L mutation that is known to confer resistance to pyrethroids ( Table 2). If all six of the mutations in the VGSC gene included in this study confer some levels of resistance to pyrethroids (kdr resistant alleles) then Southern CA populations have a higher proportion of susceptible alleles relative to central CA populations which is consistent with a previous report (13). Although Exeter 2014 and 2018 samples had different overall genetic profiles, they both contained a greater number of pyrethroid susceptible kdr alleles than the Clovis population. Presence of higher frequencies of kdr susceptible alleles may have been a factor leading to elimination of Ae. aegypti in Exeter in 2014, while other neighboring central CA districts struggled to manage dispersal of their Ae. aegypti populations.
These findings confirm that implementation of intensive control efforts can eliminate Ae. aegypti locally. However, the  replacement of the Exeter 2014 population with individuals related to the Southern CA population that had a higher proportion of pyrethroid resistance alleles adds to the difficulty of adulticide control that could be further exacerbated when they hybridize to CA Central Valley populations that have an even greater number of kdr alleles. Elimination strategies require intensive efforts and resources, including public consent and establishment of a strong public relations program. Anticipating challenges resulting from the dispersal of resistant populations throughout California will be critical to comprehensive and sustainable control strategies. This is the first record of southern CA-like Ae. aegypti genetic cluster found north of Kern county [latitude 35.21N; (25)]. The results from our study suggest that the population of southern CA Ae. aegypti continues to expand its range northwards and will hybridize with the existing Clovis population, corroborating previous indications of this possibility presented in Lee et al. (25). Additional statewide surveillance on the geographic distribution of Ae. aegypti genetic clusters combined with socio-environmental parameters may inform the nature and potential mechanisms of Ae. aegypti dispersal pathways within CA. Our relatively low-cost SNP genotyping assay or a similar approach could be a cost-effective way to screen a larger number of samples than current whole genome sequencing approaches and would be a useful tool to elucidate genetic mixing, origin of introductions, and pyrethroid resistance status of local populations.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the NCBI BioProject repository with accession number PRJNA725510.

AUTHOR CONTRIBUTIONS
YL, GA, GL, and AC conceived experimental design. EK, LM, and GA conducted SNP genotyping analysis. AR-W, T-YC, and KK conducted genomic DNA library preparations for whole genome sequencing. YL, MC, and TC conducted genome data analysis. EK, LM, KB, CG, RR-C, KS, LC, and EB conducted sample collection or arranged sample acquisition and provided the sample metadata. All authors contributed to the article and approved the submitted version.

FUNDING
We acknowledge funding support from the Pacific Southwest Regional Center of Excellence for Vector-Borne Diseases funded by the U.S. Centers for Disease Control and Prevention (Cooperative Agreement 1U01CK000516), CDC grant NU50CK000420-04-04, the Southern IPM Center (Project S21-002) as part of USDA National Institute of Food and Agriculture Crop Protection and Pest Management Regional coordination Program (Agreement No. 2018-70006-28884), the USDA National Institute of Food and Agriculture (Hatch project 1025565), UF/IFAS Florida Medical Entomology Laboratory fellowship to T-YC, and Florida Department of Health (Contract CODQJ). The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the funding agencies.