{'fraction_non_alphanumeric': 0.03624520244337532, 'fraction_numerical': 0.0223525595978161, 'mean_word_length': 5.282730514518594, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genomic Shifts, Phenotypic Clines, and Fitness Costs Associated With Cold Tolerance in the Asian Tiger Mosquito

Abstract 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).
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; 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 highthroughput 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 . 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.

Genotype-Environment Associations
The pool-target sequencing of 221 candidate genes for nine populations (PEN, PAH, PKE, LUN, HAN, FOS, SYP, The proportion of detected putative adaptive patterns was 6.8% among candidate regions, with an average F ST = 0.774, and 0% among noncandidate regions supporting a prominent role of directional selection over genetic drift in candidate genes detected as putatively adaptive. 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 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 Genomic Shifts, Phenotypic Clines, and Fitness Costs · https://doi.org/10.1093/molbev/msac104 MBE 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 temperaturelatitudinal 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 F ST ≈ 0.4), we may underestimate the proportion of outlier SNPs with BayeScan that uses population-specific F ST 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 Sherpa et al. · https://doi.org/10.1093/molbev/msac104 MBE 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 F ST genome-scan approach in A. albopictus native range . 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.

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 N G0 . 500 individuals except for MEI (N G0 , 50), LUN (N G0 ≈ 200), and SEL (N G0 , 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).

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 . We retained uniquely aligned reads and with a minimum mapping quality of 20 using SAMTOOLS v1.10 ). 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.

Library Preparation
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 SureSelect XT Target Enrichment System for Illumina Paired-end Sequencing Library protocol vD0. Genomic DNA (200 ng) was fragmented using the SureSelect XT 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.

Variant Calling
Reads were mapped to the reference genome Aalbo_primary.1 using BWA-MEM v0.7.5 . We removed secondary alignments and duplicates using Picard tools FixMateInformation and MarkDuplicates and SAMTOOLS v1.10 ). 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 persample 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 F ST ), 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 (F ST ) 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 Genomic Shifts, Phenotypic Clines, and Fitness Costs · https://doi.org/10.1093/molbev/msac104 MBE (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 . 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 (logtransformed 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.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.

Acknowledgments
For the collection of samples, we warmly thank: Qi Wang, Qi Qi, Lin Ding, Yixi Fan, Yingjie Liu, Siyun Zhou, Pingying Teng, Yuji Wang, Xiaoqing Zhang, and Pingying Teng (Southern Medical University); Phoutmany Thammavong, Phonesavanh Luangamath, Nothasin Phommavan, Somphat Nilaxay, and Vaeky Vungkyly (Institut Pasteur du Laos). This research received funds from the National Key R&D Program of China (2020YFC120100 and 2016YFC1200500), from Kasetsart University (FF[KU]14-64), from the Laboratoire d'Ecologie Alpine LECA (Université Grenoble Alpes). We thank support from the Environmental and Systems Biology (BEeSy) structure (Université Grenoble Alpes) for laboratory rearing. We thank the support of Franck Delbart for cold exposure, which was conducted within the Long-Term Socio-Ecological Research (LTSER) platform Lautaret-Oisans, a site of the European Research Infrastructure eLTER and of AnaEE-France. S.S. was supported by postdoctoral fellowships from the Université Grenoble Alpes and CNRS. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contributions
S.S. and L.D. designed and coordinated the study, with support from J.-P.D. S.K., I.H.I., X.G., J.S., S.B., S.M., T.C., X.-G.C., and X.Z. coordinated and performed the sampling. S.S., J.T., and T.G. performed laboratory rearing and phenotypic measurements. S.S., L.D., and F.L. performed molecular biology work. S.S. analyzed the data. S.S. and L.D. wrote the manuscript. All authors have read, contributed substantially, and approved the final version of the manuscript.

Data Availability
The raw sequences (fastq) analyzed in the present study are available at the European Nucleotide Archive repository (http://www.ebi.ac.uk/ena) under study accession number ERP136066, and individual accession numbers are provided in supplementary table S2, Supplementary Material online. The ddRADseq SNP dataset (1,901 ddRADseq SNPs), allele frequencies at candidate genes (exon-capture data), phenotypic data, and wing morphometric data in the form of wing shape and size variables, are available at Dryad Digital Repository (DOI: https:// doi.org/10.5061/dryad.rxwdbrvbg).