Genomic Shifts, Phenotypic Clines, and Fitness Costs Associated With Cold Tolerance in the Asian Tiger Mosquito.

Climatic variation is a key driver of genetic differentiation and phenotypic traits evolution, and local adaptation to temperature is expected in widespread species. We investigated phenotypic and genomic changes in the native range of the Asian tiger mosquito, Aedes albopictus. We first refine the phylogeographic structure based on genome-wide regions (1,901 double-digest restriction-site associated DNA single nucleotide polymophisms [ddRAD SNPs]) from 41 populations. We then explore the patterns of cold adaptation using phenotypic traits measured in common garden (wing size and cold tolerance) and genotype-temperature associations at targeted candidate regions (51,706 exon-capture SNPs) from nine populations. We confirm the existence of three evolutionary lineages including clades A (Malaysia, Thailand, Cambodia, and Laos), B (China and Okinawa), and C (South Korea and Japan). We identified temperature-associated differentiation in 15 out of 221 candidate regions but none in ddRAD regions, supporting the role of directional selection in detected genes. These include genes involved in lipid metabolism and a circadian clock gene. Most outlier SNPs are differently fixed between clades A and C, whereas clade B has an intermediate pattern. Females are larger at higher latitude yet produce no more eggs, which might favor the storage of energetic reserves in colder climate. Nondiapausing eggs from temperate populations survive better to cold exposure than those from tropical populations, suggesting they are protected from freezing damages but this cold tolerance has a fitness cost in terms of egg viability. Altogether, our results provide strong evidence for the thermal adaptation of A. albopictus across its wide temperature range.

Introduction

Temperature is one of the main determinants of species distribution as it impacts physiological functions and ultimately population growth (Pörtner et al. 2006). The tropical regions show low seasonal variation and high ambient temperatures compared with higher latitudes marked by strong seasonality, daily thermal fluctuations, and cool temperatures. These climatic characteristics impact the thermal tolerance of ectotherms (Sunday et al. 2011) and species having a latitudinal distribution are expected to present fitness shifts due to local thermal adaptations between climatic regions (i.e., ecotypes), or more gradual changes along the climatic gradient (i.e., phenotypic clines). The main fitness-related traits associated with climatic gradients reported in insects include body size (temperature–size response) (Forster et al. 2012; Horne et al. 2018) and survival at low temperature: diapause and cold tolerance (Denlinger 1991; Andersen et al. 2015; Sinclair et al. 2015).

Diapause is one of the crucial ecological adaptations explaining the success of insects to colonize high latitudes and expand their range. However, sudden temperature changes often precede diapause induction by photoperiod, and frost events can occur in spring or end of summer in temperate regions. Thermal stress is a potentially important selective agent and climatic uncertainty could favor the evolution of cold tolerance independently from the diapause programming. The relationship between diapause and cold tolerance has been a subject of controversy (Denlinger 1991) and both distinct and convergent gene expression profiles were reported in diapausing and nondiapausing mosquito eggs (Poelchau et al. 2013). The tiger mosquito, Aedes albopictus, is an invasive species originating from Asia, where populations occur across a wide latitudinal gradient from Malaysia to Japan (fig. 1). Although A. albopictus populations from temperate regions possess a photoperiodically inducible diapause of eggs, tropical populations do not (Hawley 1988; Lounibos et al. 2003), and another invasive mosquito species (A. aegypti) that lacks diapause was thought to be restricted to tropical regions. However, a cold-tolerant A. aegypti temperate population exhibiting photoperiod-induced egg-dormancy was recently described in Argentina (Fischer et al. 2019). Egg survival after cold exposure depends on the origin (higher for temperate populations), the temperature and the duration of exposure, and both cold acclimation and diapause increase cold tolerance (Hanson and Craig 1994; Thomas et al. 2012; Tippelt et al. 2020). Nonetheless, the diapause in A. albopictus is not an immediate reaction to environmental cues but is hormonally activated by photoperiod stimulus received by the female that produces diapausing eggs (Denlinger and Armbruster 2014). Recent experiments showing similar survival rates after short cold exposure between diapausing and nondiapausing eggs from the same population suggest that cold tolerance can express independently from the diapause programming (Tippelt et al. 2020).

Fig. 1.

Distribution of sampling localities. Sampling point colors represent the country of origin. In addition to the 41 newly sampled localities, we added seven localities (in italic) from a previous survey (Sherpa, Blum, Després 2019). All sample sites were genotyped at ddRADseq loci. For phenotypic data and exon-capture data, the color of the star indicates: both (black), only phenotypic (violet), and only exon-capture (red) data. Sampling characteristics in supplementary table S1, Supplementary Material online.

Here, we examined the temperature–size and cold tolerance responses in A. albopictus native range in East Asia using common garden experiment. In order to characterize the pattern of variation (i.e., clines vs. ecotypes), we sampled nine populations representative of the temperature gradient. We specifically addressed the impact of unpredictable and unfavorable temperature change in colder environments by comparing the cold tolerance of populations reared in nondiapause-inducing conditions for both temperature and photoperiod. Body size was estimated by wing size, a reliable proxy in mosquitoes (Van Handel and Day 1989), and cold tolerance by egg viability after short cold exposure. Because the mechanisms leading to increased cold tolerance may have a negative effect on reproductive success, we measured two components of fitness in control conditions (fecundity and egg viability), in order to evaluate the cost of cold tolerance when rearing conditions simulate an optimal season for growth and reproduction.

In order to determine to what extent phenotypic variation corresponds to the evolutionary history of populations, we first refined the phylogeographic structure of A. albopictus in East Asia (Kotsakiozi et al. 2017; Sherpa, Blum, Després 2019) using 1,901 double-digest restriction-site associated DNA sequencing (ddRADseq) single nucleotide polymophisms (SNPs) and large geographical sampling (191 individuals from 41 populations; fig. 1, supplementary table S1, Supplementary Material online). Studies using high-throughput sequencing (HTS) revealed a genetic toolkit underpinning diapause and cold tolerance in A. albopictus. Transcriptome sequencing (RNAseq) quantified the patterns of gene expression under diapause and nondiapause-inducing conditions (Poelchau et al. 2013; Huang et al. 2015) and landscape genomics analyses based on reduced-representation HTS (ddRADseq) identified outlier SNPs between populations from low, medium, and high latitudes (Sherpa, Blum, Després 2019). We searched for genomic signatures of selection using 51,706 SNPs located in 221 candidate genes involved in sensory perception, stress resistance, diapause, and/or cold tolerance, using pool-target DNA capture from nine populations distributed across the mean temperature gradient in Asia.

Results

Population Genetic Structure

The 125 newly sampled A. albopictus yielded 228 million ddRADseq reads with 1.8 ± 0.8 million reads per individual on average (supplementary tables S1 and S2, Supplementary Material online). Ten individuals with a high amount of missing data (≥60%) were removed from the dataset. Seventy-six previously analyzed individuals from Malaysia, China, and Japan were included (Sherpa, Blum, Després 2019). The reference-based variant calling included filtering on linkage and low frequency variants that could affect population genetic structure analyses. The final dataset contained 1,901 SNPs for 191 individuals from 41 populations (fig. 1), with average coverage 37 ± 16× and percentage of missing data 17 ± 8% (supplementary table S2, Supplementary Material online).

Three distinct genetic groups were identified by all population genetics analyses (fig. 2). The first genetic cluster includes individuals from all southern localities in Laos, Cambodia, Thailand, and Malaysia (clade A; tropical climate), the second those from China and Okinawa (clade B; subtropical to temperate climate), and the third all the individuals from Japan (except those from Okinawa) and South Korea (clade C; temperate climate). Interestingly, most localities in Northern China showed 30–50% admixture with the northernmost clade C but were all clustered together within clade B in the MCMC tree. Population SYP from Okinawa Island (Japan) was genetically and climatically closer to Chinese populations and belongs to clade B.

Fig. 2.

Population structure of Aedes albopictus in Asia. A: PCA showing genetic variation among individuals on the two first PCs. Individual points are colored according to the country of origin. B: Geographic distribution of the three genetic clusters inferred from the ADMIXTURE analysis. Pie charts represent individual ancestry coefficients (K = 3) averaged by sampling sites. C: coancestry matrix between each pair of individuals inferred using fineRADstructure, including the fineRADstructure MCMC tree with posterior probabilities (pp) >0.90 shown, and ADMIXTURE ancestry coefficients of all individuals at K = 3. The coancestry matrix represents individual (upper diagonal) and average (lower diagonal) coancestry coefficients. Individual font colors as in PCA. Boxes highlight the three clades.

Phenotype–Environment Associations

Phenotypic traits were obtained for nine A. albopictus populations distributed along the temperature gradient in Asia, reared during two generations (G) in control conditions and using five replicates (fig. 1, supplementary figs. S1 and S2, Supplementary Material online). Wing geometric morphometrics were performed for 314 individuals including 219 newly digitalized pictures from the nine populations (SEL, FOS, MEI, HAN, PKE, LUN, SYP, HKM, NGT) and 95 previously analyzed from Malaysia, China, and Japan (PBA, PSU, BAY, HAZ, PAY, SAI) revealed weak wing shape variation among populations. The PCA on Procrustes shape coordinates differentiated five populations from Japan and China on PC1 (25.51%), and populations from South China and Malaysia on PC2 (14.76%) (supplementary fig. S3, Supplementary Material online). Wing size was significantly different between populations and clades: individuals of clade C were on average larger (1.27 ± 0.00) than those of clade A (1.18 ± 0.01) and B (1.22 ± 0.00), but the factor with the strongest effect on size variation was mean annual temperature (MTEMP) (table 1, supplementary figs. S4 and S5, Supplementary Material online). Wing size was negatively correlated with MTEMP, when considering only the 219 individuals of the present study (Pearson’s R = −0.63, t = −12.10, df = 217, P < 0.001), the 95 individuals previously analyzed (Pearson’s R = −0.71, t = −9.44, df = 93, P < 0.001), or all individuals (Pearson’s R = −0.63, t = −14.35, df = 312, P < 0.001) (fig. 3).

Table 1.

Results of ANOVA Testing the Effect of Mean Temperature, Clade, and Population (nested inside clade) Factors on Phenotypic Variance.

Response	Explanatory	Df	Sum of squares	Mean square	F-value	P-value
Fecundity	Mean temperature	1	73.3	73.3	0.2368	0.6300743
	Clade	2	17342.9	8671.5	28.0210	1.368e−07***
	Clade/Population	5	10572.0	2114.4	6.8325	0.0002316***
	Residuals	30	9283.9	309.5
Viability	Mean temperature	1	0.55799	0.55799	100.7983	6.03e−11***
	Clade	2	0.02214	0.01107	1.9996	0.1536265
	Clade/population	4	0.18781	0.04695	8.4818	0.0001167***
	Residuals	29	0.16054	0.00554
Cold tolerance	Mean temperature	1	1.47710	1.47710	131.3819	2.737e−12***
	Clade	2	0.19398	0.09699	8.6270	0.001148**
	Clade/population	6	0.12134	0.02022	1.7988	0.134328
	Residuals	29	0.32604	0.01124
Wing size	Mean temperature	1	0.258734	0.258734	280.8585	<2.2e−16***
	Clade	2	0.034583	0.017291	18.7700	2.093e−08***
	Clade/population	11	0.082164	0.007469	8.1082	2.107e−12***
	Residuals	299	0.275446	0.000921

** P < 0.01, *** P < 0.001.

Fig. 3.

Linear regressions between morphometric, phenotypic, and climatic variables. (A) Relationship between wing size and MTEMP, including previously analyzed, newly analyzed, and all specimens. (B) Relationship between G1 fecundity (mean number of eggs per female) and MTEMP. (C) Relationship between G2 viability (egg hatching rate) in control conditions and after short cold exposure (cold tolerance) and MTEMP. (D) Relationship between mean fitness (G1 fecundity × G2 viability) and cold tolerance. Point colors according to genetic groups: A (tropical): green; B (subtropical to temperate): orange; C (temperate): blue.

Fecundity (mean number of eggs per G1 female) was not significantly correlated with MTEMP (Pearson’s R = 0.04, t = 0.27, df = 37, P = 0.789) (table 1 and fig. 3). Populations of clade B had a significantly higher mean fecundity (68.7 ± 7.8) than populations of clade A (36.2 ± 7.0) and C (23.8 ± 2.5), with the highest for the two Chinese populations HAN (99.1 ± 7.7) and FOS (87.9 ± 6.9) (table 1, supplementary figs. S4 and S5, Supplementary Material online). Egg viability (proportion of G2 eggs reaching fourth instar larva) was not different between clades A (80 ± 4%) and B (72 ± 3%) but significantly lower for populations of clade C (53 ± 2%) in control conditions. After cold exposure (−4 ± 1 °C), eggs of clade C survived much better (65 ± 5%) than those of clades B (24 ± 2%) and A (13 ± 3%) (table 1, supplementary figs. S4 and S5, Supplementary Material online). Mean temperature had a larger effect on egg viability than “Population” and “Clade” factors (table 1). Egg viability was positively correlated with mean temperature (Pearson’s R = 0.78, t = 7.26, df = 35, P < 0.001) but this relationship became negative when eggs were exposed to cold (Pearson’s R = −0.84, t = −9.23, df = 37, P < 0.001) (fig. 3). Eggs from temperate climate had similar survival whether exposed or not to cold, whereas egg viability approached 100% for tropical populations in control conditions but was strongly reduced by short cold exposure (fig. 3, supplementary figs. S4 and S5, Supplementary Material online). It results in nonlinear negative relationship between cold tolerance and mean fitness (G1 fecundity × G2 viability in control conditions; five replicates per population): populations have a good fitness, a good cold tolerance, or low but never high values for both traits (fig. 3).

Genotype–Environment Associations

The pool-target sequencing of 221 candidate genes for nine populations (PEN, PAH, PKE, LUN, HAN, FOS, SYP, NGT, HKM) yielded around 283 million reads, with 31.4 ± 3.7 million reads per sample on average (fig. 1, supplementary tables S3 and S4, Supplementary Material online). Among 268,720 variant positions (50,600 ± 6,000 SNPs per sample on average), we retained only SNPs with no missing data, resulting in 51,706 SNPs with a mean coverage of 388× per sample (supplementary table S4, Supplementary Material online). We detected signatures of selection among candidate regions using two genome-scan methods (BayeScan, latent factor mixed model [LFMM]) and the strength of correlation between allele frequency and MTEMP (correlation coefficient R > 0.75 or R < −0.75). BayeScan and LFMM account for the confounding effect of neutral differentiation among populations. However, to further evaluate whether putative adaptive patterns can result from genetic drift across the nine populations, we performed the same analyses with noncandidate genomic regions (ddRAD loci) in the same nine populations (subset of the initial ddRAD dataset: 548 ddRAD loci, 4,844 SNPs, supplementary table S6, Supplementary Material online).

Among candidate genomic regions, BayeScan identified 311 SNPs with a Q-value < 0.0001, among which 89 showed a significant correlation between allele frequency and MTEMP (fig. 4, supplementary table S5, Supplementary Material online). LFMM identified 2,800 SNPs with a Q-value ≤ 0.0001, of which 105 SNPs were also detected by BayeScan, including 76 SNPs showing a significant correlation between allele frequency and MTEMP (supplementary fig. S6 and table S5, Supplementary Material online). Among noncandidate genomic regions, LFMM identified 120 SNPs with a Q-value ≤ 0.0001, of which 25 showed a significant correlation between allele frequency and MTEMP, but none of the SNPs were identified by BayeScan (minimum Q-value = 0.0005; supplementary table S6, Supplementary Material online). The mean differentiation (average BayeScan locus-specific FST overall nine populations) was FST = 0.432 among candidate regions and FST = 0.406 among noncandidate regions. The proportion of detected putative adaptive patterns was 6.8% among candidate regions, with an average FST = 0.774, and 0% among noncandidate regions supporting a prominent role of directional selection over genetic drift in candidate genes detected as putatively adaptive.

Fig. 4.

Genomic signatures of selection. Variation in allele frequencies among the nine populations from clades A, B, and C for the 76 outliers SNPs detected by both BayeScan and LFMM using a FDR of 0.0001 and significantly correlated to MTEMP (correlation coefficient ≤−0.75 or ≥0.75). Predicted effect at each position: S, synonymous mutation; NS, nonsynonymous mutation; I, intron; 5UTR: 5′ untranslated transcribed region. For the 15 genes: gene ID, gene description based on Aalbo_primary.1 annotations, gene function based on UniProt, and source of candidate genes. The list of 221 genes can be found in supplementary table S3, Supplementary Material online.

The 76 outlier SNPs were distributed across 15 candidate genes in coding (55 in exons) and noncoding regions (4 in 5′-UTR and 17 in introns) (fig. 4). We detected 12 nonsynonymous SNPs in five genes. Most outlier SNPs were differently fixed between clades A and C. Clade B had an intermediate pattern: more similar to clade A for four genes and to clade C for 11 genes, but with strong variation in allele frequencies between the two mainland populations HAN and FOS and the SYP island population. The putatively adaptive genes comprise seven genes involved in metabolic processes including three genes encoding fatty acid synthases, three genes involved in cellular processes, five genes involved in signal transduction including the circadian clock gene encoding the developmental protein timeless and two genes involved in sensory perception (fig. 4).

Discussion

The invasion dynamics of A. albopictus has been largely investigated in the past decade (Goubert et al. 2016), but genetic variation in its native Asian range received less attention. Recent work using next-generation sequencing (ddRADseq) resolved both large and fine-scale population genetic structure, and suggested differential thermal adaptation between populations located in tropical (Malaysia), subtropical oceanic (Southern China), and temperate (Japan) climates corresponding to three distinct genetic groups (Kotsakiozi et al. 2017; Sherpa, Blum, Després 2019; Sherpa, Blum, Capblancq et al. 2019). In the present study, we emphasize the advantages of coupling (1) a thorough population and genetic sampling to better describe global phylogeographic patterns, (2) fitness measurements in common garden experiments, and (3) genotype-by-environment association focusing on top candidate genes to better characterize the shape of adaptive patterns (ecotypes vs. gradual clines). Based on a representative sampling of 41 populations distributed along the latitudinal gradient in Asia, we confirm the existence of three evolutionary lineages. Allele frequencies at putatively adaptive SNPs from nine populations support three thermal ecotypes. The analysis of phenotypic variation further shows clinical variation with increasing wing size and cold tolerance along the temperature gradient. We describe genomic and phenotypic variation along the gradient of MTEMPs, used as a proxy of the latitudinal temperature-photoperiod selective pressure acting on the female to produce eggs resistant to winter conditions. Minimal temperature during the coldest month or the coldest quarter (BIO6 and BIO11), which determine the survival of eggs, are highly correlated with mean temperature (correlation coefficient >0.98, supplementary fig. S7, Supplementary Material online).

In accordance to Bergmann’s rule, individuals are larger in temperate than in tropical environments. This pattern has been found in several insect species (e.g., Bai et al. 2016; Rohner et al. 2018; Wilson et al. 2019; Wonglersak et al. 2020) including A. albopictus (Sherpa, Blum, Després 2019) but the reverse pattern was also reported, so that no general pattern of size response to temperature has yet emerged in insects (Shelomi 2012). Most studies on insect size along temperature gradients use natural populations or museum specimens: the observed pattern of larger size with cooler temperature could be related to extended development time in cold environments. In the present study, the temperature–size response is not related to developmental temperature because all populations were reared in the same control conditions. Larger size might provide A. albopictus with an evolutionary advantage in colder climate by favoring the storage of energetic reserves, allowing the females to starve for longer periods of inactivity during adverse conditions. In line with this hypothesis, caloric protein and lipid contents were shown to correlate with body size (wing length) in A. albopictus, and reserves at emergence determined median adult survival times, which ranged from 16 days at 32 °C to 100 days at 17 °C (Briegel and Timmerman 2001). Interestingly, three out of our top putatively adaptive genes are involved in lipid biosynthesis processes, further supporting the hypothesis of different energy storage strategies selected across the temperature range studied. Furthermore, more reserves may allow large females to invest cryoprotectant compounds into eggs. Indeed, the survival of eggs from the cold-adapted lineage (clade C) was not affected by short cold exposure (similar survival for eggs exposed or not), whereas this treatment strongly negatively impacted the survival of eggs from the two other clades. This suggests that eggs from the temperate clade C are more protected from freezing damages. Differential accumulation of cryoprotectant compounds in the eggs of temperate and tropical populations has been demonstrated in other insects, for example, carbohydrates in the locust (Wang et al. 2010) and triglycerides in the mosquito A. aegypti (Mensch et al. 2021). Further analyses of the metabolites that could act as cryoprotectants in the eggs of temperate populations of A. albopictus are needed to better understand the mechanisms preventing freezing damages. Whatever the mechanism involved, our study shows that cold adaptation is costly, as cold-tolerant eggs have much lower survival rate in control conditions.

Relevant fitness-related phenotypic traits include reproductive success and survival. Here we characterized the mean fitness of the population as fecundity × viability in control conditions. Despite previous analyses demonstrating that wing length is an accurate predictor of fecundity in A. albopictus (Armbruster and Hutchinson 2002), we observed no relationship between female size and fecundity in the present study (supplementary fig. S8, Supplementary Material online). O’Donnell and Armbruster (2009) found larger size-specific fecundity in Malaysia than in Japan, but concluded that the observed variation could be due to either local adaptation or genetic drift (strong population-level differentiation). Accordingly, we also found the tendency for higher size-specific fecundity in Malaysia compared with Japan (supplementary fig. S8, Supplementary Material online), but with large population effects suggesting that fecundity is determined by intrinsic population factors (supplementary figs. S4 and S5, Supplementary Material online). It could include population genetic diversity (He) as Japan and Malaysia harbor lower He than China (Sherpa, Blum, Capblancq et al. 2019) but the estimation of genetic diversity indices from our SNP dataset is difficult because it has been generated from lab-reared individuals that may have impacted their relatedness.

Local thermal adaptation was assessed using a measure of cold tolerance, that is, the survival of eggs at low exposure temperature, and by searching for signatures of selection in a large set of candidate genes. We found a higher survival rate of individuals at low exposure temperature when they originate from cold compared with tropical environments as found in other studies on A. albopictus (Hanson and Craig 1994; Thomas et al. 2012; Tippelt et al. 2020). The viability in control conditions is positively correlated with MTEMP, suggesting that mean population fitness could be negatively impacted by the ability to survive to short exposure at low temperature during the optimal season for reproduction and growth. Individuals from temperate regions, which exhibit larger body size, do not have higher fecundity, but have higher cold tolerance (i.e., local thermal adaptation). Because individuals were reared in noninducing diapause conditions, our results suggest an independent response from diapause programming, although the same mechanisms (e.g., crypoprotectants accumulation) might be involved in diapausing conditions. Accordingly, a recent study showed no difference in survival after cold exposure between diapausing and nondiapausing A. albopictus eggs (Tippelt et al. 2020).

Among the targeted candidate genes, we retained 15 genes containing outlier SNPs that show both strong population differentiation and correlation with mean temperature. None of the noncandidate regions (ddRAD loci) were identified as outliers using the same false discovery rate (FDR) threshold. This supports the role of selection rather than drift alone to drive the strong divergence observed across populations along the temperature-latitudinal gradient at these 15 genes. Allelic frequencies at outlier SNPs change at two main temperature-latitudinal points: between populations of clade C and clades A/B (∼16 °C) and between populations of clades C/B and A (∼23 °C) supporting the presence of three distinct thermal ecotypes (fig. 4, supplementary fig. S6, Supplementary Material online). Because the average genetic differentiation among the nine populations is high (BayeScan FST ≈ 0.4), we may underestimate the proportion of outlier SNPs with BayeScan that uses population-specific FST to control for population demographic history (i.e., genetic drift). We thus highlighted genes exhibiting extreme allelic frequency changes along the temperature gradient. Clade B has a particular position both for allele frequency and for phenotypic clinal variation (intermediate between clades A and C). Such pattern is expected when there is local adaptation in the presence of continuous gene flow, with selection acting mainly on major genes directly involved in cold resistance (generating fixed differences), whereas adaptive but costly phenotypes, that are determined by multiple genes (i.e., quantitative traits) and potentially under complex trade-off, such as size or fecundity, produce more clinical profiles.

Outliers include seven genes involved in metabolic processes, comprising the synthesis of fatty acids, which could reflect the changes in metabolism associated with diapause and cold hardiness mechanisms. We also identified the circadian clock gene encoding the developmental protein timeless involved in biological rhythms and especially in insect diapause (Denlinger and Armbruster 2014) supporting adaptation of A. albopictus across the latitudinal gradient, in which seasonality and mean temperatures affect biological processes. Differential expression patterns between diapausing and nondiapausing conditions in A. albopictus have been demonstrated for eight of the detected genes (Poelchau et al. 2013; Huang et al. 2015). We also confirm some putatively adaptive patterns previously identified by a ddRAD FST genome-scan approach in A. albopictus native range (Sherpa, Blum, Després 2019). This previous study used a FDR of 0.05 and screened genes within a 20 kb distance of outlier SNPs identified between three populations from low, mid, and high latitudes. In the present study, we sequenced the protein-coding regions of this list of genes, using a FDR of 0.0001, and nine populations to better represent the temperature-latitudinal gradient. Among the 32 genes considered as good candidates in the previous study (Sherpa, Blum, Després 2019), we confirm here that six are highly differentiated among nine populations and strongly associated with temperature.

Among the 15 temperature-associated detected genes, only five had nonsynonymous mutations significantly correlated with mean temperature (fig. 4), as expected for causal mutations. All the other outlier SNPs were either synonymous, or in noncoding regions. Indeed, although we targeted the exons of genes, we also obtained sequences from introns and some 5′-UTR regions. Interestingly, three over five genes show both putatively adaptive (nonsynonymous mutations) and differential expression patterns, and including a 5′-UTR mutation for LOC109429583. Changes in genetic regulatory networks rather than in proteins appear to be particularly important for rapid evolution driven by climate change (reviewed in Franks and Hoffmann 2012) and selection targeting regulatory elements were shown in several cold-adapted populations of insects including flies (von Heckel et al. 2016; Huang et al. 2021; Wiberg et al. 2021) and ants (Cicconardi et al. 2020). Sequencing the regulatory regions of the targeted genes may have increased our chance to detect genomic signatures of selection, since many of them show differential expression in diapause versus nondiapausing conditions in A. albopictus (Poelchau et al. 2013; Huang et al. 2015). However in most cases the regulatory regions of these candidate genes discovered by transcriptomic analyses are unknown.

Genome-scan from ddRADseq data can provide a first step to produce a list of candidate genes to be further investigated through target sequencing, and the combination of a candidate gene approach and expression profile can narrow the large set of candidate genes underlying adaptation. Assessing the adaptive significance of the individual alleles of each candidate would require the implementation of reciprocal translocation in the field (not recommended for highly invasive mosquitoes of public health concern), selection experiments (Pardo-Diaz et al. 2015), and/or functional validation using gene editing or RNAi silencing. Here we show by combining common garden experiments and genomic analyses at noncandidate and candidate regions, that divergent genotypes across the temperature-latitudinal gradient have different phenotypic traits and differential fitness at cold exposure. Indeed, although the cold exposure we performed does not fully capture the climatic conditions that may act simultaneously in the field to shape the observed phenotype, it allows a direct estimation of the fitness impact of temperature change on differently adapted genotypes.

Characterizing the variation of phenotypic and genetic traits among native populations is crucial for understanding the colonization dynamics of populations outside the native range (Sherpa and Després 2021); the larger the sampling, the higher the probability of identifying the precise source population for global scale invasion studies. Here, we do not find evidence for recent reintroductions (i.e., large scale, human-aided migration) within the native range of A. albopictus. Signatures of admixture between the three lineages suggest substantial gene flow between populations but it overall does not blur the phylogeographic pattern. Putatively adaptive genetic variation patterns support three thermal ecotypes corresponding to the three lineages. As studies inferring colonization routes suggest that the niche of invaded ranges correspond to the climatic niches of the source populations (Kotsakiozi et al. 2017; Sherpa, Guéguen, et al. 2019), we provide useful information for predicting invasive phenotypes at the global scale based on the genetic, phenotypic, and ecological traits of source populations and thus better evaluating the preadaptation hypothesis in A. albopictus.

Materials and Methods

Field Sampling, Laboratory Rearing, and Data Collection

Sampling

Specimens of A. albopictus were collected in 41 localities between 2017 and 2020 in seven countries: 6 in Malaysia, 2 in Thailand, 3 in Cambodia, 8 in Laos, 13 in China, 6 in South Korea, and 3 in Japan (fig. 1, supplementary table S1, Supplementary Material online). Two types of sampling were performed: adults for 32 localities and eggs for 9 localities. Adults were collected from the field or from laboratory colonies and kept in 75% ethanol at −20 °C or in Silicagel. Eggs were either field F0 (China, Japan), F1 reared in laboratory to obtain enough eggs (Malaysia) or to break diapause (Japan), or F11 from laboratory colonies (Laos) (supplementary table S1, Supplementary Material online).

Phenotypic Data

We considered all the eggs as the G0 generation and reared them during two generations in standard laboratory conditions (T °C: 27 °C ± 1 °C; humidity: 70% ± 5%) and long-day photoperiod treatment (14L:10D) (supplementary fig. S1, Supplementary Material online). G0 sampled sizes were NG0 > 500 individuals except for MEI (NG0 < 50), LUN (NG0 ≈ 200), and SEL (NG0 < 50). We measured three phenotypic traits: wing morphometrics (size, shape), fitness (fecundity, viability), and cold tolerance (viability to cold exposure). Fecundity was estimated by the average number of eggs laid by G1 females using five replicates of ten females (R1–R5) (supplementary fig. S1, Supplementary Material online). Viability was estimated as the average hatching rate of G2 eggs in control conditions using five replicates (R1s–R5s). Cold tolerance was estimated as the average hatching rate of G2 eggs after short cold exposure using five replicates (R1c–R5c). The exposure temperature and duration were determined on the basis of winter temperatures in the temperate region of eastern Asia and of hatching rates measured for tropical populations of A. albopictus after cold exposure (Thomas et al. 2012). Indeed, exposing eggs to a low temperature and long period might have resulted in observing nearly identical low/null hatching rates for all tropical to subtropical populations, whereas we aimed at testing the continuous relationship between cold tolerance and mean temperature. Although some eggs from tropical strains hatch after exposure at −5 °C when exposed during 4 h, exposure duration can be up to 1 day when the exposure temperature is only −2 °C (Thomas et al. 2012). We thus used an exposure temperature of −4 °C and an exposure duration of 4 h, for which we expected ∼50% of hatching rate for temperate populations and <20% for tropical populations (Thomas et al. 2012). Exposure was performed in a climatic chamber (Aralab), in a small thermobox to limit thermal inertia. The exposure temperature was constant and controlled using two humidity and temperature data logger (REED R6020) (supplementary fig. S2, Supplementary Material online). Wing morphometrics was performed for G1 females (supplementary table S1, Supplementary Material online).

Genetic Data

For population genetics purposes, we performed ddRADseq genotyping for one to four individuals for each of the 41 localities for a total of 125 individuals. For localities sampled as eggs, we genotyped four G1 females. We collected published ddRADseq data for seven other localities (76 individuals) from a previous survey in Malaysia, China, and Japan (fig. 1, supplementary table S1, Supplementary Material online; Sherpa, Blum, Després 2019). Concerning adaptive genomic patterns, we performed pool-targeted DNA sequencing (exon capture) in candidate genes for the cold-response in A. albopictus for nine sampling localities (FOS, HAN, HKM, LUN, NGT, PAH, PKE, PEN, SYP). For egg-sampled localities, we used 50 G1 adults (25 males and 25 females) when the number of G1 females was enough (i.e., not for SEL and MEI) (supplementary table S1, Supplementary Material online). For Malaysia, pool-targeted DNA sequencing was performed using 50 field-collected adults (25 males and 25 females) from two populations (PEN, PAH).

ddRAD Sequencing

Library Preparation

Total genomic DNA was extracted using the cetyl trimethyl ammonium bromide chloroform/isoamyl alcohol protocol (Doyle and Doyle 1987). We used the ddRAD protocol from Sherpa, Blum, et al. (2019). Genomic DNA (200 ng) was digested using 30 units each of SbfI-HF and MspI restriction enzymes (New England Biolabs, USA). Digestions (3 h at 37 °C) were ligated to Illumina P1 and P2 adapters by 90 cycles of digestion at 37 °C (2 min) and ligation at 16 °C (4 min) with 400 units of T4 ligase (New England Biolabs, USA), followed by a heat-inactivation step at 65 °C (10 min). Digested-ligated products were purified with 1.5:1 ratio Ampure XP beads (Beckman Coulter, USA). We selected DNA fragments 190–600 bp long with a Pippin-Prep 2% gel cassette (Sage Sciences, USA). Each library was amplified with 2 μL of size-selected DNA, 0.2 mM of dNTPs, 0.15 μM of each PCR primers, 3% dimethyl sulfoxide, and 0.4 U of Taq Phusion-HF (New England Biolabs, USA). The PCR conditions were: initial denaturation at 98 °C (10 min), 12 cycles of 98 °C (10 s), 66 °C (30 s), 72 °C (1 min), final extension period at 72 °C (10 min). Final libraries were purified with 1.8:1 ratio Ampure XP beads and sequenced on Illumina NextSeq500 (2 × 150 bp; Fasteris, Switzerland).

Variant Calling

Sequences were trimmed to equal length of 110 bp using the BBMAP package v37.33 (Bushnell 2014). Paired-end reads were mapped to the A. albopictus reference nuclear genome Aalbo_primary.1 (GenBank accession GCA_006496715.1) using BWA-MEM v0.7.5 (Li and Durbin 2009). We retained uniquely aligned reads and with a minimum mapping quality of 20 using SAMTOOLS v1.10 (Li et al. 2009). STACKS v2.0 (Rochette et al. 2019) was used for variant calling and SNP genotyping. We removed called genotypes from per-sample reads depth <5× and required a locus to be in ≥70% of the individuals. We removed SNPs with minor allele frequency <0.01 and retained one random SNP per locus.

Pool-targeted DNA Sequencing

Candidate Genes

We targeted 221 candidate genes based on AaloF1.2 annotation (supplementary table S3, Supplementary Material online). We included 32 genes showing strong population differentiation in A. albopictus native range (Sherpa, Blum, Després 2019), 121 genes involved in sensory perception (Lombardo et al. 2017), 14 genes involved in stress resistance (Rinehart et al. 2007), and 54 genes involved in the diapause of insects, including the insulin-foxo signaling pathway, circadian clock genes, and enzymes involved in carbohydrate and lipid metabolism (Poelchau et al. 2013; Denlinger and Armbruster 2014; Huang et al. 2015).

Genomic DNA of 50 individuals was extracted from using the PureGene kit (Qiagen, Germany). The capture of genomic regions of interest was performed using the SureSelect target enrichment system (Agilent Technologies). Capture probes were designed based on AaloF1 genome assembly (GenBank accession GCA_001444175.1) and consisted of 32,616 overlapping RNA probes of 120 bp. Capture was performed with the SureSelectXT Reagent kit (Agilent Technologies) following the SureSelectXT Target Enrichment System for Illumina Paired-end Sequencing Library protocol vD0. Genomic DNA (200 ng) was fragmented using the SureSelectXT Low Input Enzymatic Fragmentation kit (Agilent Technologies), ligated to adaptors, purified using Agencourt AMPure XP, and amplified by PCR using Herculase II Fusion DNA polymerase (Agilent Technologies). Libraries were hybridized to biotinylated baits and purified using Dynal MyOne streptavidin beads (Invitrogen). Captured DNA fragments were amplified, purified, and multiplexed before sequencing on an Illumina NextSeq500.

Reads were mapped to the reference genome Aalbo_primary.1 using BWA-MEM v0.7.5 (Li and Durbin 2009). We removed secondary alignments and duplicates using Picard tools FixMateInformation and MarkDuplicates and SAMTOOLS v1.10 (Li et al. 2009). Variant calling and SNP genotyping were performed using the Genome Analysis Toolkit GATK v4.1.9.0 (Van der Auwera et al. 2013). The HaplotypeCaller and SelectVariants tools were first used to produce per-sample VCF and select SNPs. After checking the distribution of VCF annotation values, we removed SNPs with QualByDepth QD < 2.0, StrandOddsRatio SOR > 3, FisherStrand FS > 60, MappingQuality MQ < 20, and coverage DP < 20, using the VariantFiltration and SelectVariants tools. The genotype of each sample for each position was obtained using the GenotypeGVCFs tool to include nonvariant sites (i.e., sample identical to reference genome). Finally, we merged per-sample VCFs and removed SNPs with missing data.

Genetic Analyses

Population Structure

Population structure was investigated using ddRADseq data. We first evaluated genetic variation among the 191 individuals using a principal component analysis (PCA) in adegenet (Jombart 2008). We estimated the most likely number of K genetic clusters using ADMIXTURE v1.3 (Alexander et al. 2009). Different values of K ranging from 1 to 10 were analyzed, and the best value was selected based on the 10-fold cross-validation errors. We used fineRADstructure v0.3.2113 (Malinsky et al. 2018) to infer population structure via shared ancestry among individuals. The coancestry matrix was constructed from the sorted 1,901 loci and individuals were assigned to populations with a burn-in period of 2,000,000 and 5,000,000 Markov chain Monte Carlo (MCMC) iterations and a thinning interval of 10,000. The MCMC tree was constructed from the default parameters.

Genomic Signatures of Selection

We searched for signatures of selection in the nine populations using the targeted exon-capture poolseq dataset (221 candidate genes). The detection of signatures of selection was first assessed using a population-based approach. We used the Bayesian algorithm implemented in BayeScan v.2.1 (Foll and Gaggiotti 2008). This method simulates an ancestral population split into diverging subpopulations evolving through genetic drift (estimated by population-specific FST), and for each locus, compares the posterior probability of a model including drift only to that of a model including both drift and selection, to estimate the probability of that locus being under selection. The MCMC algorithm was run for 1,000,000 iterations following a burn-in period of 100,000 with the proposal distributions for parameters adjusted by 20 short pilot runs of 5,000 iterations. The final sample size was 100,000 iterations (one sampled every 10). We then used a genotype–environment associations (GEA) method testing association between allele frequency of each of the 51,706 SNPs and climatic variation. GEA were performed using the LFMM (Frichot et al. 2013). The latent factors estimated by LFMM are included in the model when performing GEA tests to correct for unobserved confounding variables (population structure). Based on population genetics results, we used three latent factors to correct for the independent population demographic histories of the three clades. LFMM models were run with 10 repetitions and 10,000 iterations with a burn-in period of 5,000 iterations. The P-values were adjusted by using a genomic inflation factor and the FDR was controlled using a threshold of Q-value ≤0.0001.

Our aim was to discover outliers among SNPs located in candidate genes involved in the cold or diapause response of A. albopictus. We used the MTEMP as climatic predictor, representing the temperature gradient in Asia as most BIOCLIM temperature variables were highly correlated (supplementary fig. S7, Supplementary Material online). We collected MTEMP from the WorldClim 2 database (Fick and Hijmans 2017) at 30 s arc resolution and averaged for the period 1970–2000. Although LFMM identifies correlation between genetic and environmental variation, BayeScan only uses between-population differentiation (FST) at each single locus to identify genomic regions more differentiated than the average differentiation across SNPs. We retained SNPs discovered by both methods and with allele frequency significantly correlated to MTEMP (correlation coefficient R > 0.75 or R < −0.75) as candidates for A. albopictus cold adaptation. We used SnpEff (Cingolani et al. 2012) and the VectorBase database (https://vectorbase.org) to annotate and predict the effect (amino acid change) of candidate SNPs. Gene function was assessed using Gene Ontology annotations of the Universal Protein Knowledgebase (UniProt, http://www.uniprot.org).

The two methods used to detect genomic signatures of selection take into account the confounding effect of drift on genetic differentiation across populations. However, given the strong genetic differences between lineages, we further evaluated the possible confounding effect of genetic drift by performing the same analyses using ddRAD data. We had both exons capture data (N = 50) and ddRADseq data (N = 4) for eight of the nine populations (FOS, HAN, HKM, SYP, NGT, LUN, PKE, PEN) (supplementary table S1, Supplementary Material online). Because PAH and SEL, for which we also had two individuals, were genetically, geographically and environmentally close, we combined these four individuals for the ninth population. Among all ddRADseq loci present in ≥70% of 191 individuals from 41 populations (overall filter), 548 polymorphic loci were present in ≥75% of individuals in each of the nine populations (within-population filter), with 4,844 SNPs.

Phenotypic Analyses

Wing geometric morphometrics was assessed using G1 females stored in 75% EtOH at −20 °C. The left wing was mounted on a glass microscopic slide, and photographed using an Olympus DP70 digital camera (3.1 MP) connected to an Olympus SZX12 stereomicroscope (DF PLAPO 1.2xPF2 objective). Twenty landmarks (LMs) located at vein intersections and termini were digitalized using tpsUtil v.1.76 (Rohlf 2006) and tpsDig2 v.2.31 (Rohlf 2008). Variation due to scale, orientation, and position was removed by applying a Procrustes superimposition to LMs coordinates using the Integrated Morphometric Package (Sheets 2014). We analyzed a total of 314 individuals, of which 219 were newly digitalized and 95 were from a previous study (Sherpa, Blum, Després 2019). We applied a PCA on Procrustes coordinates to evaluate wing shape variation.

In order to test for gradual changes along the climatic gradient (phenotypic clines), shifts between climatic regions (ecotypes) and intrinsic population factors, we performed ANOVA testing the effect of “Mean temperature,” “Clade,” and “Population” (nested within clade) factors on phenotypic variation for wing size (log-transformed wing centroid size), fecundity, viability, and cold tolerance. The significance of differences in mean phenotypic traits between populations and between clades was assessed using pairwise multiple comparisons tests (Tukey post hoc tests) and the relationship between phenotypic traits and mean temperature was assessed using Pearson’s correlations.

Click here for additional data file.