Spatial, Temporal, and Genetic Invasion Dynamics of Aedes albopictus (Diptera: Culicidae) in Illinois

Abstract

The spread of the Asian tiger mosquito, Aedes albopictus Skuse, throughout the United States has implications for the transmission potential of vector-borne diseases. We used a 30-yr data set of occurrence records in Illinois and developed a hierarchical Bayesian model to shed light on the patterns and processes involved in the introduction and expansion along the northern edge of the geographic range of this species. We also collected specimens from 10 locations and sequenced a segment of their mitochondrial COI genes to assess possible introduction sources and geographic patterns in genetic variation present within contemporary populations. We documented an increase in the number of observations throughout the southern and central parts of Illinois over the study period. The process through which this spread occurred is likely only partially due to local dispersal. The probability of successfully overwintering was likewise low, but both these parameters increased over the study period. This suggests that the presence of Ae. albopictus has been largely due to repeated introductions, but that in recent years populations may have become established and are leading to an increase in locally driven dispersal. There was considerable genetic diversity among populations in Illinois, with 13 distinct haplotypes present in 10 sampling locations, several of which matched haplotypes previously found to be present in locations such as Texas or Japan. Further research is needed to understand how the combination of continued propagule pressure and establishment of populations are driving the increase and expansion of this invasive mosquito along its northern distribution limit.

Biological invasions have consequential repercussions on ecosystems and the population processes of native species embedded in them. The rate at which such invasions occur and their impact on the environment is inextricably linked to anthropogenic activities, notably the increasing activity in international trade (affecting propagule pressure) and habitat degradation (affecting receptivity and potentially invasive potential). Whether a given species is likely to itself become invasive can depend on a wide variety of intrinsic factors, including its abundance in, and breadth of, the native range, closeness of contact with humans and human activities, as well as its life history and ability to respond to natural selection via its genetic architecture (Lee 2002, Hufbauer et al. 2012, Prinzing et al. 2002).

Understanding the processes that lead to successful invasions is a critical part of developing management or mitigation methods. As a general framework, we can distinguish between distinct phases in the invasion process: introduction to new areas, establishment (i.e., a positive rate of population growth in the novel environment), and further geographic spread and ecological impact following establishment (Mack et al. 2000, Shea and Chesson 2002). Typically, invasive alien species are introduced unsuccessfully many times before becoming established and can then remain established at low or unnoticeable levels for long periods before suddenly and drastically increasing in abundance and spreading (Crooks and Soulé 1999, Sax and Brown 2000, Allendorf and Lundquist 2003). Such patterns likely result from complex interactions between propagule pressure (i.e., the numbers of introduced specimens per time unit, likely eggs in the case of Ae. albopictus), environmental and demographic stochasticity, landscape heterogeneity, ecological interactions with predators, competitors, and pathogens, local adaptation, and local environmental change over time. Disentangling these factors and their impacts on lags in range expansion and population growth, and understanding to what extent adaptation to novel environments is driving range expansions, remains a challenging problem in invasion biology (Crooks and Soulé 1999).

In the case of invasions by vectors, such as mosquitoes, the ramifications of their introduction and spread extend to the realm of public health. Such introduced species can alter the transmission level of pathogens in their new range through various direct and indirect mechanisms (e.g., different levels of vector competence or human-biting preference than native mosquitoes, or via alteration of mosquito species assemblages) (Lounibos 2002, Juliano and Lounibos 2005).

A prime example of a highly successful invasive species of public-health concern is Aedes albopictus Skuse, the Asian tiger mosquito (Lounibos 2002). This species has been shown to be a highly competent vector under laboratory conditions for a wide variety of arboviruses, including dengue, Zika, chikungunya, West Nile, and La Crosse encephalitis viruses (Talbalaghi et al. 2010). While there are a number of dengue epidemics where Ae. albopictus was the putative primary vector, including in Gabon, Indian Ocean islands, as well as in southern Europe (Gasperi et al. 2012, Bonizzoni et al. 2013), it has been noted that those outbreaks lack the large numbers of dengue hemorrhagic fever cases associated with epidemics caused by the more anthropophilic Ae. aegypti L. (Diptera: Culicidae) (Gratz 2004, Lambrechts et al. 2010). In the case of chikungunya virus (CHIKV), Ae. albopictus has become a major concern and vector in many parts of the world (Bonizzoni et al. 2013). This illustrates the public health threat associated with an aggressive vector that can attain high densities and spread to novel environments, where its presence has the capacity to change the selective forces acting on endemic or emerging viruses.

Aedes albopictus is thought to have originated as a sylvatic, zoophilic mosquito in South-East Asia, and thrives in rural and peridomestic habitats in tropical and temperate zones in this region (Mogi et al. 2015). Over the past half century this species has spread and become established over a large part of the globe (Bonizzoni et al. 2013), in large part facilitated by the international trade in used tires (Reiter and Sprenger 1987, Benedict et al. 2007). From the 1980s onwards, it was introduced to and spread through large parts of the United States, following its introduction in Texas, and also in Europe, where it occurs throughout the Mediterranean region and is likely expanding northward (Benedict et al. 2007, Scholte and Schaffner 2007, Kraemer et al. 2015c). Early on, spread appears to have been closely related to the presence of interstate highways, supporting the notion that human-mediated dispersal may have been a main driver of range expansion (Moore and Mitchell 1997).

There are a number of traits that have been associated with invasiveness in mosquitoes. For instance, species with drought-resistant eggs that can survive transportation over distances and times likely have higher probabilities of being introduced to new environments (increasing propagule pressure), though not necessarily of spreading or becoming invasive in those new environments. An ability to thrive in human-dominated environments, however, appears to be related to invasion success (Juliano and Lounibos 2005). This may be because an association with such disturbed habitats reflects an ability to occupy an otherwise relatively empty niche (e.g., small containers in peridomestic environments).

In an analysis of the presence of Ae. albopictus in the United States between the period of 1995–2016, 40 states had reported occurrence records, with the densest distribution of the counties in the southern and eastern states up to a line of about 40°N (Hahn et al. 2016).

An open question, with potential repercussions for vector-borne disease transmission patterns, is to what extent the range of Ae. albopictus will continue to expand northward, and through which (combination of) genetic, ecological, and environmental processes (e.g., climate change and increasing urbanization) (Kraemer et al. 2019).

Early studies suggested that Ae. albopictus populations could overwinter in areas with mean January temperatures warmer than 0°C (Nawrocki and Hawley 1987). Areas where daily mean January temperatures are between 0 and −5°C would be amenable to late summer expansions only. Isothermal lines at 0° and −5° January temperature could potentially then delineate the northern range expectations of Ae. albopictus (Armstrong et al. 2017), The climatic conditions also depend on the geographic origin of the Ae. albopictus population, as evidenced by considerable differences in temperature tolerances within the native region (Nawrocki and Hawley 1987). It is also possible that populations can evolve greater cold hardiness and continue to adapt to local conditions (Hanson and Craig 1995, Medley et al. 2019). The latter is of interest, because habitat suitability models based on the native range of this species failed to predict a significant proportion of the occurrence records in North America related to the northward and westward expansion (Medley 2010), suggesting that ecological niche of this species has shifted, possibly via adaptation, as a result of factors related to the invasion process. For instance, introduced populations may have undergone a bottleneck and experienced founder effects. Alternatively, a continuous influx of new introductions, likely mediated by human transport (Medley et al. 2015), of different backgrounds may provide genetic material for adaptation to occur.

To develop a better understanding of the establishment and spread (and the potential for future spread) of this invader along its current northward range limit, we first use the historic record to define where and when Ae. albopictus invaded the state of Illinois and then use statistical analyses and genetic evidence to shed light on the processes that have influenced its historic invasion and establishment in the state.

Methods

We first assembled county-level data from a combination of published literature (Kraemer et al. 2015a,b; Hahn et al. 2016, 2017), the Illinois Natural History Survey Insect Collection, and data reported to the Illinois Department of Public Health (J. Swanson, unpublished) to indicate when Ae. albopictus was observed and reported for 102 counties in the state of Illinois, from the first year of observation in 1986 to 2015 (Fig. 1). Then, we developed data as indicators of the potential for 1) introduction of Ae. albopictus into a county, 2) conditions that allow for overwintering, 3) diffusion from adjacent counties, and 4) observation opportunities. We used these data in a multivariate hierarchical Bayesian statistical model to determine the processes that best explained the pattern of annual, county-level observation of Ae. albopictus. In addition, we analyzed a sequence of mitochondrial DNA (mtDNA) from collections in 10 locations and assessed the genetic structure of Illinois Ae. albopictus (Fig. 2). For spatial data processing and map development, we used ArcGIS v. 10 (ESRI, Redlands, CA). For statistical analysis and model implementation, we used the No-U-Turn sampler (NUTS) (Matthew and Gelman 2014), an adaptive form of Hamiltonian Monte Carlo sampling (Neal 2011, Marwala 2012), in Stan (Carpenter et al. 2017). Convergence was assessed using Gelman-Rubin diagnostics (Gelman and Rubin 1992), with values less than 1.2 suggesting convergence.

Occurrence Records

We developed a county-level dataset of observations of Ae. albopictus by year and by county with records from multiple sources (Smith 2019). An observation was defined as any record of any number of any life stage of Ae. albopictus collected in Illinois where the year and county of collection were both available. The collection years under consideration were from 1986 to 2015. All counties used approximately the same adult trapping protocol, as provided by the Illinois Department of Public Health. Potential sources included state and local health departments, CDC Mosquito Net, Illinois mosquito abatement districts, the Illinois Natural History Survey Medical Entomology Laboratory, published peer-reviewed papers, and individual scientists. All observations were compiled first in a list, with each row being an observation. Then, the list was reconfigured into a matrix made up of the total number of observations summarized as the 102 Illinois counties by the 30 possible years.

Estimate of Overwintering Potential

The potential for Ae. albopictus to overwinter was estimated as a combination of the average January temperature and the size of the human population in a particular county. Two overwintering potential interval values were developed and then added together to form an index of overwintering potential that provided an estimate of the degree to which winter temperatures might affect overwintering.

Data for the 30-year average temperature in January was downloaded from the CHELSA climate project (Fig. 2). The CHELSA data are at 30 arc second resolution, and the ‘normals’ are for the period from 1979 to 2013 (Karger et al. 2017). We used the ArcGIS tool to summarize by zone to calculate the mean value from the gridded PRISM data of the mean January temperature for each county. Then, an ordinal measure was developed based on temperature above or below 0°C, whereby higher values represent the places where overwintering is more likely, based on Nawrocki and Hawley (1987). A value of ‘0’ was given to the 11 counties where mean January temperature was colder than −5°C and there was little chance of winter survival; ‘1’ to 48 counties between −2°C and −5°C; ‘2’ to 23 counties between 0°C and −2°C; and ‘3’ to the warmest 20 counties between 1.7°C and 0°C.

Due to the heterogeneous distribution of the human population of Illinois, there is a high likelihood of urban survival of mosquitoes due to the urban heat island effect and the increase in potential overwintering habitats, such as peri-urban tire dumps and catch basins. Therefore, the population size of each county was determined from the U.S. Census Bureau count of human population in 2000. The Illinois county population distribution is strongly skewed by the over five million residents in Cook County (Chicago, IL), which is five times larger than the next largest county, DuPage. Of 102 counties, about 80% have a population of fewer than 100,000 people. Based on this distribution, and the idea that larger urban areas are more likely to harbor favorable microclimates and habitats during cold winters (Cardo et al. 2014, LaDeau et al. 2015), an ordinal variable was developed whereby: a value of ‘0’ was given to 84 counties with population less than 100,000; ‘1’ to nine counties with 100,000 to 200,000; ‘2’ to eight counties with 200,000 to 1,000,000; and ‘3’ to Cook County, alone.

Spatial Diffusion and Detection

To account for the potential for diffusion from a neighboring county that had a prior observation of Ae. albopictus, we used the U.S. Census Bureau county adjacency file (U.S. Census Bureau 2018) and incorporated a term into the statistical analysis that measured ‘1’ if any adjacent county had previously observed Ae. albopictus and ‘0’ otherwise. To account for the variable likelihood of there being resources and interest to collect Ae. albopictus in a county, we used a measure based on the number of mosquito traps that were set out for West Nile virus surveillance. As a measure of collection effort, we used the average number of locations where mosquito traps were set out per year during the years from 2004 to 2013 as reported to the Illinois Department of Public Health and based on the West Nile Virus surveillance protocol recommended by that entity. An ordinal variable was developed whereby a value of ‘0’ was given to seven counties where there were no locations with mosquito traps set for mosquito surveillance; ‘1’ to 44 counties with from more than zero to two trap locations; ‘2’ to 30 counties with from 2 to 10 trap locations; and ‘3’ to 21 counties with from 10 to 240 trap locations. It is assumed that the detection index will also capture the relative effort to identify adult mosquitoes collected in the traps, which we believe to correspond closely with the number of trap locations.

Statistical Model

The first stage of the analysis involved the development of time-series maps of the observations of Ae. albopictus to review general patterns. Through the review of spatial correlations and data distribution of the indicator variables, as well as an understanding of Ae. albopictus biology, we selected the factors we expected to be most influential for determining the probability of presence, namely the overwintering and detection indices and diffusion potential (presence in neighboring county), for inclusion in the Hierarchical Bayesian model. Second, we formulated the model (Fig. 3) which estimates the probability P_i,t of Ae. albopictus presence in a given county (i) and a given year (t). Our parameters of interests are therefore P^O (the probability of overwintering), P^D (probability of diffusion), and P^E (the probability of detection given presence). The first level (data/observation level) in our hierarchical model is [Y| P^O (w^O), P^D, P^E(w^E)]. For succinctness, we suppress the dependence of P on w. In the second level (process level), since P^O and P^E are already modeled as a function of the calculated indices (X^O and X^E) with unknown parameters w^O and w^E, this level is given by [P^D|α,β]. We close the hierarchy by giving the priors of [w^O], [w^E] [α,β]. The model assumes that the probability of presence is a function of the independent probabilities of overwintering $Pi,tO$ and diffusion (P)$Pi,tD$⁠,

and then that observations of Ae. albopictus are Bernoulli random variables with probability $Pi,t∗Pi,tE$⁠, where E indicates the detection effort index of county i. Both the overwintering and effort functions take a logistic function; for J representing either overwintering as O or effort as E, given the calculated indices $Xi,tJ$⁠, we assume a logistic function for $Pi,tJ$⁠:

The joint posterior distribution of all unknown model parameters is then defined as

where [A|B] represents the conditional distribution of A given B, [A] represents the distribution of A, and α and β are the shape and scale parameters for a Beta distribution representing diffusion probability from a neighboring county when at least one neighboring county has previously observed Ae. albopictus. The explicit formulations of each term in the posterior distribution are

where Y_k indicates the observation of Ae. albopictus in a neighboring county k. The prior distributions of the w^J are assumed to be Gaussian, with mean 0 and standard deviation 100 (Table 1). The prior distributions of α and β are assumed to be mixtures of U ~ uniform (0,1) and V ~ gamma (1,20), where α = U * V and β = (1 − U) * V.

Mosquito Collections, Sequencing, and Population Genetic Analysis

Adult female Ae. albopictus were collected in the summer of 2017 in 10 locations across central and southern Illinois, using BG-Sentinel or gravid traps (Table 2). For the sites where BG-Sentinel traps were used, collections were performed on a biweekly basis between July 25 and October 12, with three traps per city. Dry ice was added to each trap to act as an additional cue. Collected specimens were transported to the laboratory in a cool box containing dry ice and were subsequently identified to species based on morphological characteristics (Darsie and Ward 2005). Samples from other locations were provided by public health and mosquito abatement districts and were obtained via routine West Nile virus surveillance programs using gravid traps. The choice of sites was related to an Ae. albopictus surveillance effort in counties in Illinois without previous occurrence records (Stone et al. in preparation) and by submissions of specimens from public health partners that resulted in a sufficient number of females for the genetic analysis.

DNA was isolated from a leg of each mosquito using the Phire Tissue Direct PCR Master Mix kit (Thermo Fisher Scientific, Waltham, MA). Obtained DNA was PCR-amplified using a MyFi Mix (Bioline USA, Inc., Taunton, MA) and two sets of primers flanking the mitochondrial cytochrome c oxidase subunit 1 (CO1) gene: 1454F (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and 2160R (5′-TAA ACT TCT GGA TGA CCA AAA AAT CA-3′); 2027F (5′-CCC GTA TTA GCC GGA GCT AT-3′) and 2886R (5′-ATG GGG AAA GAA GGA GTT CG-3′) (Zhong et al. 2013). The amplifications were performed using Applied Biosystems Veriti thermal cyclers. The PCR cycles consisted of a 94°C denaturing step for 3 min, followed by 35 cycles consisting of 30 s at 94°C, 30 s at 52°C, and 1 min at 72°C, with an extension of 6 min at 72°C. The amplicons from the two PCRs were sequenced separately using primers 2160R and 2027F on Applied Biosystems 3730xl DNA Analyzers at the Roy J. Carver Biotechnology Center at the University of Illinois at Urbana-Champaign.

After examining the electropherogram of each sequence using Unipro UGENE v1.30, a sequence from 2160R primer and a sequence from 2027F primer derived from an individual Ae. albopictus were joined at the overlapping region between reads to obtain a 1321 nucleotide-long sequence. The resulting gene sequences were aligned using Clustal Omega (Sievers et al. 2011) and then edited manually using AliView (Larsson 2014). Unique sequences were resequenced for confirmation. The number of variant sites, haplotype diversity (h), and nucleotide diversity (π) were estimated for each of the different locations using DnaSP version 6.10.04 (Rozas et al. 2017). As a measure of genetic distance between sampling sites we calculated pairwise F_ST values using Arlequin (Excoffier and Lischer 2010). The distm function of the geosphere package in R was used to calculate the shortest geographic distances between sets of coordinates following the Vincenty (ellipsoid) method (Hijmans 2016). A Mantel test performed in Arlequin was then used to test for a potential correlation between genetic and geographic distances.

To visualize genetic relationships between samples, and potentially aid in making inferences about the biogeographical history of populations, we created haplotype networks using POPART (Leigh and Bryant 2015). The haplotype network was inferred using a median-joining network method (Bandelt et al. 1999). We used a Bayesian computational method to investigate potential genetic structuring of Ae. albopictus individuals into clusters with limited gene flow. Specifically, we used the software BAPS, which implements a stochastic partition-based mixture model (Corander et al. 2008). We used the option of clustering with linked loci to perform the genetic mixture analysis for haploid sequence data and provided the sampling population as additional information for each individual. Singletons, haplotypes that occurred only in a single specimen, were removed for this analysis. We performed 20 individual runs. The program determines the optimal number of genetic partitions or clusters, k, for each run and merges these based on their likelihoods.

Results

Spatial and Temporal Dynamics

A total of 2,387 records of Ae. albopictus in the state of Illinois were identified after removing all potential duplicate records. The first report of Ae. albopictus was in 1986 in St. Clair and Jefferson counties. Ae. albopictus was also observed in Cook County in 1987, in the northeastern part of the state (Moore et al. 1988), but observations since that time have been limited to 2016 and the years since (Fig. 1). By 2015, 40 of 102 counties had reported observations of Ae. albopictus at some point in time. The years 2000, 2011, and 2014 each had eight or nine counties with observations, making them the highest years in this regard. Most of the occurrences (2,337, 97.9%) were in the southern two-thirds of the state (between about 37° to 41° N latitude).

Three counties had consistent observations over time during the study period compared to the other counties, where observations were intermittent. Peoria County, at the northern edge of the range of counties with observations, saw Ae. albopictus consistently every year starting in 1997, due to consistent ovitrap surveillance by the regional office of the Illinois Department of Public Health. St. Clair County, in the southwest part of the state was one of the two counties (with Jefferson) that were among the first observations in 1986, and it had 12 more subsequent years with observations. Macon County, in the center of the state, had 14 yr of consistent observations starting in 2002. The other 45 counties had from 1 to 3 yr with observations.

The posterior of the overall probability of Ae. albopictus presence,$P(Pi,tO,Pi,tD,Pi,tE|Yi,t)$⁠, is shown in Fig. 4. There is a high likelihood of presence in Central Illinois starting in 2001 and continuing through 2015, while Southern Illinois has a higher likelihood of presence in the decade of 1991–2000. Considering the individual components of this probability, the probability of overwintering was estimated to remain low throughout the state and the probability of diffusion, though increasing over time, was also low. The model also predicted that a moderate-to-high level of effort in detection was likely to find Ae. albopictus, if present, but counties with very low detection effort had a low probability of detection.

Population Genetics

After alignment, we obtained sequences of 1,321 nucleotides of the CO1 mitochondrial gene for 245 individual Ae. albopictus females. Within this region, for the individuals we analyzed, there were 11 polymorphic sites, of which nine were parsimony informative, resulting in a total of 13 haplotypes (Supp Table S1 [online only]). The overall haplotype (h) and nucleotide diversity (π) among all sampling locations were 0.694 (±0.023) and 0.0015 (±0.00005), respectively. Among sampling locations, haplotype diversity ranged from 0.197 in specimens collected from Marshall to 0.754 for individuals from Decatur. Nucleotide diversity ranged from 0.00039 in Marshall to 0.00163 in Decatur (Table 2). There was no apparent association between geographical locations and genetic diversity among samples.

The haplotypes present in these Illinois populations can be divided into two major groups, consisting of two common haplotypes and a number of smaller ones (Fig. 5). These two major haplotypes (Hap_1 and Hap_2) are found in almost all locations (Hap_1 not being present in Albion being the exception), though the proportional abundance of each varied dramatically. For instance, haplotype 1 was more abundant in the western parts of the state, while haplotype 2 was predominant along the eastern border of the state. This is reflected in the outcomes from the clustering analysis of samples, which resulted in nine partitions associated with the larger haplotypes (Fig. 6). Among these (and excluding haplotypes present in only single individuals), all haplotypes clustered independently, save for Hap_8 and Hap_10, which were joined in a single cluster (no. 1). Cluster 5 (Hap_3) is prevalent in the central and southern locations, while cluster 9 (Hap_6) can be found in the northern locations. Two of the clusters were unique to specific locations: cluster 2 (Hap_12) to Mt. Carmel and cluster 3 (Hap_13) to Cass.

Certain of the haplotypes present in these communities correspond to those found by Zhong et al. (Zhong et al. 2013), which putatively sheds insight on the potential sources of these populations. For instance, our haplotype 1 corresponds to their h37, which they found in populations from Italy, New Jersey, and Texas. Our Hap_2 corresponds to their h25, a haplotype they found in a population from Japan. Our Hap_3 is equivalent to their h55, found in a population from Texas, while Hap_6 corresponds to h39, present in Italy and Texas. A haplotype represented by a single specimen, Hap_7, corresponded to h17 (found in Taiwan, Italy, California, Texas, and Hawai’i). All other haplotypes we found had no complete match to any of the haplotypes identified by Zhong et al. or any other perfect matches according to GenBank BLAST searches.

We tested for isolation-by-distance, i.e., whether genetic differentiation between locations was related to their geographic distance. A comparison of pairwise F_st (a measure of population structure) and geographic distance between sites suggests no such correlation exists among these populations (Mantel test: r = 0.105, P = 0.245).

Discussion

Our investigation of occurrence records and the contemporary genetic diversity of Ae. albopictus documents the geographic spread and increase of this species in Illinois over a 30-yr period, and highlights the diversity of genetic backgrounds that has resulted from or contributed to this range expansion. We find overall that while the probability of introduced populations successfully overwintering was low, the probability with which presence in a given county was estimated to be due to dispersion from adjacent counties (i.e., a result of establishment) increased over time. The results also highlight the low probability of detecting this species throughout this period. The pattern of observations over this period shifted from a number of seemingly isolated observations in the southern part of the state to more frequent, widespread, and consistent observations in the central areas of the state. This change in the pattern of observations raises several important questions regarding the invasion of Illinois by Ae. albopictus. One question is whether the shift in observations from the south to the center of the state reflects range expansion and spread, or a shift in the actual presence of this species within Illinois. In other words, are the absences of observations in the south through most of the second half of the study period reflective of true absences? The notion that this species has spread and increased its presence in Illinois is also evident from the posterior probabilities of presence (Fig. 4), which largely follows the occurrence records. These estimates likewise suggest a low probability of presence in the final two time periods, even though the model is accounting for low detection efforts in these areas. One plausible possibility is that basing the detection effort and probability on current West Nile virus surveillance capacity fails to capture the true detection rate for Ae. albopictus and its variation over time. Inclusion of more recent observations that fell outside of the study period used for this modeling exercise may provide a different picture for these regions. For instance, in collections by the INHS Medical Entomology Lab in 2016 and 2017 in 18 south-eastern counties in Illinois, Ae. albopictus was found to be present throughout (Kim & Stone, unpublished). Thus, the low or intermittent level of capacity to collect Ae. albopictus, in part of the state suggests that the lack of observations in later years may provide relatively little information regarding the continuing presence of this species. In that case, the cumulative observations taken over the entire 30-yr period may provide a better overview of the actual presence.

The probability of overwintering in Illinois was found to be fairly low, even for counties below the 0° isotherm of Nawrocki and Hawley (1987). However, there is likely to be confounding related to this finding. The counties in Southern Illinois with the highest overwintering index values also tended to have low values for the detection effort index, resulting in an observation bias. Hence, there is a likelihood that overwintering potential was high, but observations were not made. In contrast, some areas of Central Illinois, which were above the 0° isotherm, saw repeated observation of Ae. albopictus, likely due to the high levels of detection effort of local public health departments in those areas. This highlights the importance of homogeneity in data collection, which was not possible for this study due to its retrospective nature.

There was a low, but increasing, probability of diffusion of Ae. albopictus between counties. These mosquitoes do not independently travel far as adults (Maciel-de-Freitas et al. 2006) and instead rely on anthropogenic modes of transport (Medlock et al. 2006, Beebe et al. 2013, Eritja et al. 2017). If such anthropogenic transport from county to county is not a major source of spread, other potential modes of introduction (e.g., tires, plants) should still be prioritized in biosecurity measures. It is notable, though, that the parameter estimates for the probability of diffusion, while low overall, rise during the study period. This, together with the overall increase in counties with observations during the study period, suggests that some of these populations may have become established and spread via shorter range dispersal modes, rather than solely being the result of an increase in propagule pressure. Both short-range dispersal and habitat suitability for Aedes spp., however, are responsive to heterogeneity in the landscape (Russell et al. 2005, Maciel-de-Freitas et al. 2010, Murdock et al. 2017). With a finer-scale and more extensive data set on Ae. albopictus presence and abundance, the effects of such local fine-scale variation on establishment and spread could be investigated further.

A pattern of repeated introductions of Ae. albopictus into the state of Illinois appears to be quite distinct from others that have been described, for instance, in southern France, where a clear dispersal wave was observed and with ‘jumps’ (introductions in areas not adjacent to areas that had previously been established) rarely being successful (Roche et al. 2015). The statistical model presented here did not consider the probability of introduction, as we found that there were no data capable of producing a realistic introduction index. Locations with facilities for the importation, collection, storage, and processing of used tires are of particular interest in assessing the potential for the transport of Ae. albopictus and for the possibility of enhancing the ability of the mosquitoes to breed locally (Baumgartner 1988, Yee et al. 2010). Therefore, we had considered the use of the Illinois database on tire pile remediation to indicate counties in which there was a high probability of introduction. However, we found that this was negatively correlated with the detection of Ae. albopictus. This is likely due to the fact that the database only records tire piles that have been removed by the state, not those that are actively importing tires or that have tires remaining. Road networks, indicated by the presence of an interstate highway, were also not associated with observations. This is likely because the presence of the mosquito in surrounding states is not known. Other posited introduction sources, such as lucky bamboo imports (Madon et al. 2002), do not provide readily available data for analysis. Thus, we chose to focus on the movement and overwintering within the state.

We found a surprisingly high level of genetic variation in mtDNA and diversity of haplotypes within Illinois, possibly a consequence of the relatively large section of the COI gene we sequenced here (Zhong et al. 2013, Goubert et al. 2016), although in their assessment of genetic diversity in samples from across the globe Zhong et al. (2013) identified 66 haplotypes in total, with individual sites harboring between 2 and 11, compared to 13 here. While there was apparent geographic separation between the haplotypes, with different genetic backgrounds dominating in different parts of the state, most haplotypes were present throughout with a diversity of haplotypes or clusters in each location (Figs. 5 and 6). The majority of specimens fell within two major haplotypes, though there were two additional haplotypes present at intermediate numbers, and a number of haplotypes present in very small numbers. This diversity, and the lack of an isolation-by-distance effect, suggests that multiple introductions of Ae. albopictus populations have occurred or are continuing to occur. As this analysis was based on mtDNA sequencing results for contemporary populations sampled in 2017, we are unable to disentangle the history or order in which introductions and possible establishment occurred over time, or whether displacement of any clusters is going on by newer, possibly better-adapted haplotypes. For instance, of the larger haplotypes (2, 3, and 6) all have previously been identified in collections from Texas as well (Zhong et al. 2013), where Ae. albopictus was initially introduced into the continental United States, likely from Japan (Manni et al. 2017). Although speculative, it is possible then that one or more of these haplotypes were introduced to Illinois as a result of the initial spread of Ae. albopictus throughout the central and eastern U.S. Haplotype 2 on the other hand, which is most abundant in the eastern part of the state, has previously only been found in specimens from Japan (Zhong et al. 2013). Further work (e.g., multi-year collections) will be required to understand whether this haplotype is in the process of spreading and potentially supplanting other haplotypes. The remainder of the haplotypes, save one which corresponded to one found previously in collections from California, Texas, and Hawai’i, all did not correspond to any haplotypes in the literature or GenBank. Whether these represent independent introductions, or are the result of genetic changes in individuals following establishment in Illinois is likewise unknown. More detailed analyses making use of genome-wide SNPs could possibly shed light on questions related to drift, whether selection is occurring, and to what extent admixture between populations is occurring.

It remains to be determined whether and to what extent this genetic diversity within Ae. albopictus is of epidemiological concern, beyond the possible implications for range expansion and population increases. Recent work suggests that the probability of becoming infected can vary strongly among different populations for viruses such as Zika and dengue (Paupy et al. 2001, Chouin-Carnero et al. 2016). Similar studies into the vector competence of different haplotypes to viruses of more local concern, such as La Crosse or West Nile, could provide useful insights. Understanding whether other traits that determine vectorial capacity (e.g., survivorship, host preferences) differ between these populations will be equally or more important.

Finally, given the documented range expansion and likely establishment of Ae. albopictus throughout the central part of the state, a pressing concern remains to what extent these populations will continue their northward expansion. A question that is of particular relevance to Illinois is whether Ae. albopictus will (re-)gain a foothold in the Chicago area, given the density of the human population that resides there and the density of international travel that passes through the area. To a large extent, over the course of the next decades, such dynamics will likely be driven by changes in climate and landscape (e.g., increased urbanization) (Rochlin et al. 2013, Ogden et al. 2014). Even in the absence of such environmental change, or genetic changes in mosquito populations, the question whether this species has currently reached the limits of its range, or whether this edge is still being colonized and expanded upon year to year is not entirely clear. A recent study suggests that Ae. albopictus may not yet have saturated its current ecological niche in the United States, supporting the notion of continued expansion in the coming decades (Kraemer et al. 2019). The concerns and the uncertainty regarding the extent of further range expansion by this important vector, as well as the inherent uncertainty associated with the absence of collections in areas with relatively low surveillance capacity suggest that an increase in systematic surveillance or surveys for the presence of Ae. albopictus along its northern edge should be considered, to allow for amelioration of any possible implications for vector-borne disease dynamics in the affected areas. Given the diversity of genetic backgrounds and possibility for hybridization between these populations, and of further adaptation to local conditions, incorporating routine assessments of genetic backgrounds can both aid in identifying the sources of these populations and provide more detailed insight into their possible roles in maintaining or altering arbovirus transmission cycles in recently invaded areas.

Supplementary Data

Supplementary data are available at Journal of Medical Entomology online.

Fig. S1. The posterior probabilities of overwintering estimated as a function of the overwintering index for each county. The count (white dot) indicates the value of the index, a sum of mean January temperature and county population size.

Acknowledgments

We are grateful to M. Blackshear, D. Lin, D. Patel, R. Back, K. Gochanour, S. Yates, and N. Paoletti for assistance with the research, and to the St. Clair and Jackson public health districts and Macon County mosquito abatement district for providing additional samples for the genetic analysis. This work was supported by the State of Illinois Used Tire Management and Emergency Public Health funds to C.M.S. C.H.K. acknowledges support from the Illinois Department of Public Health. Points of view or opinions expressed in this document are those of the authors and do not necessarily represent the official position or policies of the Illinois Department of Public Health. This publication was also supported by Cooperative Agreement U01 CK000505, funded by the Centers for Disease Control and Prevention to R.L.S. and M.O.R. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the Centers of Disease Control and Prevention or the Department of Health and Human Services. This work was supported by National Science Foundation award 1830312 to B.L. and R.L.S.

References Cited