The effects of genetic drift and genomic selection on differentiation and local adaptation of the introduced populations of Aedes albopictus in southern Russia

Sample collection

Mosquito samples were collected between September 17 and 19, 2019 at six locations (Fig. 1) in the Krasnodar Krai (Russian Federation). Three samples (30 mosquitoes in 1_Sochi, 17 in 2_Sochi, and 30 in 3_Sochi) were collected in Sochi, where the climate is humid subtropical without a marked winter. Three samples (42 mosquitoes in 4_Krasnodar, 25 in 5_Krasnodar, and 40 in 6_Krasnodar) were collected in Krasnodar with a more continental climate and below zero winter temperatures. Locations within a city were separated by at least 10 km, while the distance between the two cities is over 100 km. Collection sites were located in city cemeteries (excluding 3_Sochi, which was located at a rabbit farm at the forest edge). Visitors to the cemetery bring many memorial vases, which serve as reservoirs, where immature stages of Ae. albopictus develop. Sampling was performed from 6 to 7 pm, when female mosquitoes attack most actively. Adult female mosquitoes were captured using a mosquito aspirator when they landed and attempted to feed on collectors. After capturing, mosquitoes were identified using the Gutsevich, Monchadsky & Stackelberg (1970) Field Guide. Insects were stored in cold 96% ethanol.

DNA extraction and high-throughput sequencing

Mosquitoes from each collection site were pooled in a 1.5 ml tube containing 180 μl of PBS. Homogenization was performed using TissueLyser II (QIAGEN, Hilden, Germany) at 30 shakes/s for 4 min. DNA was extracted from homogenized samples using the DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany) in accordance with the manufacturer’s protocol. Short (150 bp) paired-end whole-genome (WGS) reads were obtained by Macrogen Inc. (Seoul, South Korea) using Illumina NovaSeq 6000. Average coverages were 58.4× for Sochi (1.9× per individual), 64.4× for 2_Sochi (3.8× per individual), 57.2× for 3_Sochi (1.5× per individual), 61.4× for 4_Krasnodar (1.5× per individual), 64.3× for 5_Krasnodar (2.6× per individual), and 58.4× for 6_Krasnodar (1.6× per individual). The mean GC-content was 40.28%. Read quality was checked using the FastQC software (Andrews, 2010); adaptors and low quality (phred < 20) reads were trimmed using Trimmomatic v. 0.38 (Bolger, Lohse & Usadel, 2014). Raw sequence reads for pool-seq samples are accessible by Sequence Read Archive accession numbers SRR13570421 to SRR13570426.

Read mapping, population genomics, and sweep detection

Reads were aligned in pairs to the Ae. albopictus reference genome (Refseq assembly AalbF2 accession GCF_001876365.2) using bbmap (Bushnell, 2014). Popoolation (Kofler et al., 2011) software was used to estimate population genetics parameters, such as π, $\theta$, and Tajima’s D, in pooled samples using a 1 kb sliding window, and Popoolation2 (Kofler, Pandey & Schlotterer, 2011) to estimate Fst and conduct Fisher’s exact test to compare differences in the allele frequency spectrum (AFS) between samples (1 kb sliding window, step-size = 1). Genomic regions for which significant (FDR < 0.05, Bonferroni correction) differences in AFS were found between the populations from different cities and not found within Sochi or within Krasnodar were selected.

Single nucleotide polymorphism and indel calling was performed using GATK HaplotypeCaller (McKenna et al., 2010). The variance-covariance matrix (Ω) for cluster analysis using the Bayesian hierarchical model was build using the neutral core model (Coop et al., 2010) in the BAYPASS software (Gautier, 2015). Effective population size was estimated using the poolne_estimate script (Gautier et al., 2013).

Selective sweep detection analysis was performed using SweeD (Pavlidis et al., 2013), grid = 1,000, for each pooled sample. For scaffolds over 50 Mbp grid = 10000 was used. Scaffolds shorter than 100 kb were excluded from the analysis. Genomic regions affected by sweep were selected based on the Composite likelihood ratio (CLR) greater than a 0.5% quantile in each sample. Gene was considered to be affected by sweep if it was located in a sweep region or if a sweep was located in a putative promotor and transcription binding site region (1 kb upstream). BLAST2GO (Conesa et al., 2005) was used for the functional annotation of genes located in sweeps.

Phylogenetic analysis

Filtered reads were mapped on the reference mitochondrial genome (NC_006817) separately for each sample using bbmap (Bushnell, 2014). These mapped reads were used to reconstruct haplotype sequences and frequencies in each sample using RegressHaplo pipeline (Leviyang et al., 2017) (max_num_haplotypes = sample size, sig = 0.01, num_trials = 1,000). Mitochondrial genome haplotype sequences of Sochi and Krasnodar samples were aligned using MAFFT (L-INS-i method, Katoh et al., 2005) with other available in GenBank Ae. albopictus mitochondrial genomes (Zé-Zé et al., 2020; Battaglia et al., 2016; Zhang et al., 2016). We used PartitionFinder (Lanfear et al., 2012) to identify a partitioned model of sequence evolution for mitochondrial genome alignment and Mr.Bayes (Ronquist et al., 2012) for Bayesian phylogenetic inference (MCMC ngen = 5,000,000). We checked the convergence of the runs by the program Tracer 1.7 (Rambaut et al., 2018).

Sequencing and read mapping

Whole genome sequencing produced 142.9–161.2 Gb of raw read sequence data, 97% of which were retained after quality trimming. About 71–77% of reads were successfully mapped on the Ae. albopictus genome. It resulted in 57.2–64.4× coverage per a pooled sample and 1.5–3.8× coverage per an individual mosquito within a pool. For mitogenome, the coverage was 417–479× per sample.

Differentiation between populations

Firstly, the suitability and representativity of the obtained sequence data were checked by comparing the sample size and the effective population size. Local populations 1_Sochi and 4_Krasnodar had a three-fold lower effective population size than its actual sample size (9.7 ± 1.7 against 30 individuals and 19.3 ± 3.5 against 42 individuals, respectively), probably due to the isolation that led to inbreeding. Population 3_Sochi showed a less drastic deviation between the effective population and sample sizes (20.3 ± 7.5 vs. 30). Other local populations had effective population size close to the actual sample size, according to the whole-genome allele frequency estimation using poolna_estimate.

Analysis of differentiation between the populations revealed 44.25 million polymorphic sites across the Ae. albopictus genome. More than half of them (22.59 million sites) were significantly (FDR < 0.05, Bonferroni correction) different between at least one pair of studied populations. Analysis of fixation indices (Fst) between samples revealed high level of population differentiation (Table 1). Differentiation between the three samples from Sochi was higher than between 1_Sochi and all Krasnodar samples, as well as between 2_Sochi and 4_Krasnodar or 3_Sochi and 4_Krasnodar. Based on the mean Fst values across the genome, only 5_Krasnodar and 6_Krasnodar appeared to be weakly differentiated; the observed differentiation between all other pairs of samples pointed to the high numbers of specific alleles in each population.

Cluster analysis also confirmed strong differentiation between the local populations except for 5_Krasnodar and 6_Krasnodar (Fig. 2). Despite Krasnodar local populations were separated by a similar geographical distance from each other, 4_Krasnodar mosquitoes were genetically as distant from other Krasnodar populations as from 2_Sochi mosquitoes.

Genomic regions with distinct allele spectrum

To identify genes associated with cold adaptation, we selected genomic regions, which showed significant differentiation between samples from Sochi and Krasnodar, but not between the samples collected within each of the two cities. These loci were filtered from outliers by the Fst value according to the Rosner’s test. As a result, we obtained 46 Ae. albopictus genomic regions (Table S1) with allele frequencies differing between the Sochi and Krasnodar, but not within the Sochi or Krasnodar populations.

Selective sweeps

Composite likelihood ratio distribution across the scaffolds with candidate selective sweeps for each local population are presented in Fig. 3. About 30–40% of these regions did not contain genes within the 10 kb range. Genes found in the selective sweep regions are presented in the Table 2.

Three of these genes were involved in fatty acid metabolism (LOC109422645, LOC115263932, LOC115260991, GO:0006633 fatty acid biosynthetic process). Many loci located in selective sweeps were also involved in development, in particular, in neurogenesis and eye development: homeobox protein HMX3, NXF1, DCX, neuronal PAS domain-containing protein 3, PAX-6, and rdgC (GO:0007399 nervous system development). Some of them participated both in the developmental pathways and pathogen response, such as Krüppel-like factors and runt-related transcription factors.

Phylogenetic analysis and haplotype distribution

Phylogenetic analysis based on the whole mitochondrial genomes showed that despite significant differences in the allele frequency spectrum Ae. albopictus populations from the Krasnodar Krai are closely related and share a common origin (Fig. 4). Tiger mosquitoes from Russia formed a separate clade with the posterior probability of 1. Mosquitoes in this clade appeared to be related to mosquitoes from Italy, Portugal, and Greece, their mitochondrial haplotypes belonging to the haplogroup A1a2. We didn’t find any shared mitochondrial haplotype between any pair of analyzed samples. 3_Sochi and 4_Krasnodar populations revealed a single unique mitochondrial haplotype each, while other samples revealed one major (present in 73–80% of the individuals according to PredictHaplo) and two to four less common haplotypes (Table 3). The COI region in all samples of the Sochi and Krasnodar mosquitoes was identical to the one previously found in this region (MG198595, Shaikevich et al., 2018).

The haplotypes identified in the 1_Sochi, 2_Sochi and 5_Krasnodar samples formed a distinct clade. In contrast, the rare 6_Krasnodar haplotypes (Krasnodar_6_H4 and Krasnodar_6_H3) showed more basal localization on the phylogenetic tree relative to all other Sochi and Krasnodar haplotypes (Fig. 4).

Experimental design limitations and coverage effects

We have selected Ae. albopictus populations from Sochi and Krasnodar because of distinct differences in winter temperatures between these two regions. However, it should be noted that Ae. albopictus can also be found in other areas of the Krasnodar Krai, where climatic conditions are different from those in both chosen locations. In Novorossiysk and Tuapse, located in the northern part of the Black sea coast, the average daily temperature does not differ dramatically from that in Sochi, but average precipitation is lower and annual average wind speed is higher (gismeteo.ru). It may cause differences in selective pressure and patterns of local adaptation. It is also expected that there is a gene flow between these places, which can limit the efficiency of selection. There is no data confirming that Sochi was the immediate source of mosquito expansion to Krasnodar, thus selective pressure imposed by some other factors apart from low winter temperatures could have shaped the variation patterns.

To accurately detect selection sweeps it is necessary to reduce the demographic history effects, since the bottlenecks and founder effects are hardly distinguishable from selective sweeps when working with small samples (Jensen et al., 2007). Thus, the 2_Sochi (n = 17) and 5_Krasnodar (n = 25) samples may be considered too small to efficiently detect sweeps with relatively small selection coefficients. To reduce demographical effects, we had planned to find selection sweeps in the same genomic regions in all three Sochi populations compared to all three Krasnodar populations to be sure that selection search results are reliable. We found no genomic regions of this kind in the Sochi and Krasnodar populations, although certain genes under selection in different samples belong to the same pathways or share the same function.

Furthermore, the sweep detection using low coverage pool-seq data may not be accurate, because of the wrong base calling or analyzing just one chromosome of a diploid organism (Pavlidis et al., 2013). For a reliable detection of 30% allele spectrum differences between two populations both coverage and pool-size of 50 are needed (Kofler et al., 2011). Most of our samples passed this threshold, but two (2_Sochi and 5_Krasnodar) were slightly deficient in pool sizes. Of note, for a pairwise Fst analysis, the results obtained with even lower coverage data did not drastically differ from the classical estimates (Hivert et al., 2018), so we may assume that the quality of our data was sufficient for differentiation analysis.

Tiger mosquito dispersal and bottlenecks

In a previous phylogenetic analysis of Ae. albopictus mosquitoes from the Black Sea region, no significant differentiation between local populations was found using COI and ITS2 data, probably because only samples from subtropical areas near Sochi were studied (Shaikevich et al., 2018). It should be noted that COI marker lacks the resolution to discriminate Krasnodar krai Ae. albopictus populations, as only two haplotypes were found previously (Shaikevich et al., 2018). In agreement with that, no variations were found in the COI region in the current work.

Genetic differentiation study using PoolSeq revealed that Ae. albopictus populations from Russia belong to the same cluster and clade (Fig. 4), suggesting that the populations within the region diverged as a result of in situ differentiation after a single or multiple introduction(s) from the same geographical region. Similar results were reported previously for Ae. albopictus invasion regions in Italy, Africa, and South America (Kotsakiozi et al., 2017). While clear differentiation between different geographical populations in Italy and Africa was demonstrated, its level (fixation indices) was about two times lower than that observed between Sochi and Krasnodar. In contrast, Ae. albopictus mosquitoes from Greece and Albania showed low level of divergence from each other, probably due to the lower rates of this mosquito expansion into these countries (Pichler et al., 2019). However, in another study, a different scenario for Ae. albopictus invasion in Italy was proposed. According to Manni et al. (2017), Italian Ae. albopictus populations were a result of multiple introductions and admixtures with the populations from Reunion and the United States. Somewhat controversial results were obtained in the two subsequent ddRADseq studies of Italian Ae. albopictus mosquitoes—no recent gene flow between Italy and Reunion (Sherpa et al., 2018), but admixture between the Northern Italy and China mosquitoes in Central Italy was demonstrated (Sherpa et al., 2019). The Ae. albopictus mosquitoes from Russia form a separate clade on the mtDNA phylogenetic tree. This implies a single or several introductions from the same geographic location. It is plausible that unrelated mtDNA haplotypes have been independently introduced to other seaports in the Krasnodar Krai, but it is unlikely that they have spread to the continental parts of the region, since dominant mtDNA haplotypes detected in Krasnodar differed non-significantly from the ones detected in Sochi (Fig. 4). This means that tiger mosquitoes from southern Russia descended from a group of closely related females.

Rather unexpectedly, we found no mitogenome haplotypes shared between the studied samples (Table 3). This could possibly be due to a sampling bias, but, taking into account fixation indices between populations (Table 1), it may rather be an indication of the highly restricted gene flow between populations. Average Fst in this case were higher than those observed between the populations from other invasion regions such as Europe, continental Africa (Kotsakiozi et al., 2017; Pichler et al., 2019), and Brazil (Multini et al., 2019).

It is noteworthy that certain populations from different cities appeared to be more closely related to each other than the population within the same city (Table 1, Fig. 2). It can be accounted for by multiple independent introductions of mosquitoes from the coast to the inland areas, multiple introductions from other areas, and very intensive diversification of mosquito populations in Sochi. The latter case may be exemplified by the 3_Sochi and 1_Sochi populations exhibiting high level of divergence. Local geographic conditions favored diversification, since these mosquito populations were localized to the rural areas separated by foothills, Sochi center, and other water basins. Additionally, 1_Sochi population was collected closer to the Adler seaport, while 2_Sochi and 3_Sochi were collected closer to the Sochi seaport. Thus, the differences between these populations could be a result of separate introductions (although from the same general source, as the mitochondrial sequences suggest) without subsequent merge of the established populations. The 4_Krasnodar population exemplifies how populations may become isolated even within a generally endemic area. It was sampled on a graveyard surrounded by urban areas. The effective population size analysis suggested that the population was affected by inbreeding, while fixation indices indicated that it was significantly different from other populations. In contrast, in the regions with a dense network of waterbodies gene flow increases (Medley, Jenkins & Hoffman, 2015). The similarity between the 5_Krasnodar and 6_Krasnodar populations could be explained by sharing the same water body, namely a Kuban river bend with multiple piers, which may contain small artificial reservoirs, and the proximity to the A146 highway leading to Novorossiysk. At the same time, the 4_Krasnodar population was found at a certain distance from the river transport system and closer to the M4 highway, thus allowing different introduction histories for these populations. This hypothesis fits well into the observation that Sochi mosquitoes are less closely related to each other compared that those from Krasnodar (Table 1). This may be due to the fact that Sochi populations are separated by foothills and rivers or the independent introduction of mosquitoes, sharing similar mitochondrial genomes, but different allele frequencies. Previously, lower genetic differentiation was shown for Aedes mosquitoes inhabiting large permanent water storages compared to those inhabiting the temporary ones (Paupy et al., 2004). In the urban areas, different pattern and density of oviposition sites (Huber et al., 2002) and availability of human-operated transportation means were a more important differentiation factor than geographical distance and natural barriers (Medley, Jenkins & Hoffman, 2015).

The Ae. albopictus spread across the Krasnodar Krai likely included multiple parallel transmissions of mosquitoes between the cities, as judged by the relatively low levels of differentiation observed between the samples from different cities compared to within-city differentiation. Very high fixation indices and absence of shared haplotypes between samples are compatible with severe bottleneck effects and, possibly, difficulties in local climatic adaptation due to reduced selection effectiveness.

Cold adaptation

There was no uniform pattern of genomic regions that differed in the allele frequency spectrum between the Sochi and Krasnodar samples detected. Some of the genes in the selective sweep regions identified in different local populations, however, shared common functions, such as fatty acid metabolism (GO:0006633) or neural system development (GO:0007399).

Genes involved in fatty acid metabolism are also often involved in cold adaptation (cold tolerance) in insects (Kayukawa et al., 2007; Michaud & Denlinger, 2006; Yang et al., 2018). In 1_Sochi, 4_Krasnodar, and 5_Krasnodar populations, selective sweeps were found in the genomic regions containing the genes encoding short-chain specific acyl-CoA dehydrogenase, proteins participating in the elongation of very-long-chain fatty acids, and fatty acid synthase-like proteins. Previously, these genes were found to be involved in egg desiccation resistance and lipid conservation during diapause (Reynolds et al., 2012; Urbanski et al., 2010). As the egg is the only stage in the tiger mosquito development cycle that can survive sub-zero temperatures (Kreß et al., 2016), it is expected that selection would act upon genes functioning in diapausing eggs. In addition to the genes involved in fatty acid metabolism, there were other loci affected by selection that could be involved in adaptation at the diapause stage. For example, chitin deacetylase and larval/pupal cuticle protein H1C may change the cuticle composition at the critical developmental stages (Chen et al., 2010; Dunning, 2013).

Selection acting on serine/threonine-protein phosphatases rdgC and Pax-6, as well as other proteins responsible for brain development, is consistent with the previous findings which linked cold adaptation to the thermotaxis rhodopsin pathway (Parker et al., 2015; Pavlides, Pavlides & Tammariello, 2011). In Drosophila spp., genes participating in neuronal development were under positive selection in cold-resistant species (Parker et al., 2018). Neural protection during diapausing stages may be crucial for further brain development.

Although field data indicate that Ae. albopictus expansion into the Krasnodar Krai occurred very recently (Ganushkina et al., 2016; Fedorova et al., 2018, 2019b), the complexity of local populations exceeded our expectations. While we achieved sufficient genome coverage, the lack of common haplotypes among populations, relatively high genetic distances, and potentially multiple introduction routes suggest that the number of sampling locations needs to be increased. Due to the technical and logistical limitations, the present study could not cover the whole dispersal history and cold adaptation of Ae. albopictus in southern Russia. In the coastal area, there were other locations, where Ae. albopictus mosquitoes were recorded. These locations have different climatic conditions to adapt to, such as the rainfall amount and wind speed. In the case of the samples used in the present work, different local landscapes and different distances to the seaports and other spread routes could result in high divergence between the Sochi populations and could have significantly affected our ability to identify genes associated with cold adaptation. In other conditions, 2–3 samples from a geographic area could be sufficient to analyze pool-seq data (Asgharian et al., 2015; Konorov, 2018), but in this study, the sampling coverage was obviously not exhaustive, likely due to background drift effects.