A total of 25 novel mitogenomes were included in this study. Twenty-two were from wild populations collected in Europe, Asia, and the Americas (Table 1; Figure 1). Three were from the Americas (two from Virginia and one from Brazil). Nine were from Asia: three from Thailand (one from Hang Chat district, Lampang province in the North; one from Ban Rai district, Uthai Thani province in the West; one from Phato district, Chumphon province in the South), five from Los Baños, Laguna, Philippines and one from Wakayama prefecture, Japan. Ten were from Europe: two from Tirana (Albania), two from Athens (Greece), two from Cesena and two from Pavia (Northern Italy), one from Cassino (Central Italy), and one from Reggio Calabria (Southern Italy). This study also included three adult laboratory-maintained strain mosquitoes: two from the Italian Rimini strain (Bellini et al., 2007; Manni et al., 2015), established at CAA (Centro Agricoltura Ambiente “G. Nicoli,” Crevalcore, Italy) from mosquitoes collected in Rimini, Italy, and one from the Chinese Foshan strain (Center for Disease Control and Prevention of Guangdong Province; Table 1). The study did not involve protected species and specimens were not collected at sites protected by law.

Morphological keys (Rueda, 2004) and/or PCRs with species-specific primers for internal transcribed spacer regions (ITS1 and ITS2) of ribosomal DNA (rRNA; Higa et al., 2010) were used to identify the specimens. For the Philippine samples, eggs were collected using ovitraps, and the emerging adults were reared in an insectary under standard conditions of temperature (27°C), humidity (60–80%), and photoperiod (12:12 h). Samples were preserved in 80% ethanol and stored at -20°C until DNA extraction. DNA extraction was carried out using the Wizard Genomic DNA Purification Kit (Promega) following the manufacturer’s protocol.

A primer set was designed to amplify the entire A. albopictus mitogenome in four overlapping PCR fragments (protocol I). The A. albopictus Reference Sequence (NC006817) was used to design primer sets. The coding region (nps 1-14893) was amplified by two long PCR fragments whereas the control region (nps 14894-16665) was amplified by two rather short PCR fragments (Supplementary Table 1) because of its high A+T content and the presence of repeated elements requiring distinctive PCR procedures.

Coding region long PCRs (Supplementary Table 1) were carried out in 50 μl reaction mixture containing 1X GoTaq Long PCR Master Mix (Promega), 0.2 μM of each primer and 10–20 ng of DNA template using the following PCR conditions: 94°C (2 min); 30 cycles of 94°C (30 s), 59°C (30 s), 65°C (9 min); and a final extension of 72°C (10 min). Alternatively, a set of nine overlapping PCR fragments (Supplementary Table 2) covering the A. albopictus coding region was also designed (protocol II). PCRs were carried out in 25 μl reaction with a standard reaction mix containing 1X White Buffer (1.5 mM MgCl₂), 0.2 mM of each dNTP mix, 0.6 U of GoTaq G2 Polymerase (Promega), 0.2 μM of each primer and 20–30 ng of DNA template, using the following PCR conditions: 94°C (2 min); 35 cycles of 94°C (30 s), 55°C (30 s), 72°C (2 min); and a final extension of 72°C (5 min).

PCR primers used to amplify the control region in two overlapping fragments (Supplementary Table 1) were the same in both protocols. The control region PCRs were carried out using the following PCR conditions: 94°C (2 min); 35 cycles of 94°C (30 s), 54°C (30 s), 60°C (2 min); and a final extension of 60°C (5 min) for PCR #3, and 94°C (2 min); 40 cycles of 94°C (30 s), 55°C (30 s), 60°C (1 min); and a final extension of 60°C (10 min) for PCR #4. PCR products were visualized on a 1–2% agarose gel and successful amplicons were sequenced with standard dideoxy sequencing using Big Dye v3.1 Chemistry (Applied Biosystems) on 3730xl and 3130xl Genetic Analyzer (Applied Biosystems) following the manufacturer’s protocol. Sets of 28 and 29 oligonucleotides were designed to sequence the A. albopictus mtDNA coding region starting from protocol I (Supplementary Table 3) or protocol II (Supplementary Table 4), respectively. Sequences were assembled and aligned using Sequencher 5.0 (Gene Codes) comparing them with the Reference Sequence (NC006817) from a Taipei sample, Taiwan.

Given the different copy numbers of repeated elements contained in the A. albopictus control region, PCR fragments #3 and #4 (Supplementary Table 1) yielded products of different length in different mosquitoes. The two types of tandem repeats, I and II (Figure 2), were amplified with PCR #3 and PCR #4, respectively. Overall, the size of the amplified fragments ranged from ∼1,800 to ∼2,500 bp for PCR #3 and from ∼800 to ∼900 bp for PCR #4. As for one of the two mitogenomes from the Rimini strain (mitogenome #1, Figure 1; Table 1), PCR #3 and #4 yielded products of ∼2,000 bp and ∼900 bp, respectively. The product of PCR #4 was directly cloned in the pCR2.1 TOPO vector (Invitrogen) following the manufacturer’s protocol. White colonies were PCR-screened for the insert length and desired clones were sequenced bi-directionally using the M13 universal primers. The procedure described above is inefficient for longer fragments such as that of ∼2,000 bp, therefore a different strategy was developed to sequence the region amplified by PCR #3. It was re-amplified in two overlapping PCR fragments (PCRs I and II), each with one primer within the tandemly repeated elements of type I (Supplementary Table 5). Multiple amplicons were obtained for each PCR and were cloned as described above. Only the two clones containing the longest fragments deriving from each PCR were sequenced bi-directionally. Sequences for each clone were assembled and aligned using Sequencher 5.0 (Gene Codes). Two additional internal primers were then newly designed (data not shown) to confirm the sequence obtained by cloning.

The obtained complete sequence of mitogenome #1 (accession number KX383916) was used to assemble and number, with Sequencher 5.0 (Gene Codes), the other 24 novel mitogenomes (accession numbers KX383917-35, KX809761-65). We aligned the novel sequences and the two previously published mitogenomes – #17 from China (Zhang et al., 2015) and #27 from Taiwan (NC006817) – by performing a Multiple-Sequence Alignment with the Clustal algorithm (Chenna et al., 2003) implemented by Sequencher 5.0.

Most Parsimonious (MP) trees (1000 bootstrap replications) encompassing the 27 mitogenomes were built by using MEGA7 (Kumar et al., 2016), employing the Tree-Bisection-Regrafting (TBR) algorithm (Nei and Kumar, 2000). Modelgenerator v.85 indicated for our dataset HKY+I as the best-supported model according to the AIC1, AIC2, and BIC criteria. The obtained settings were selected to infer maximum likelihood (ML) and Bayesian trees for the A. albopictus dataset. The ML tree was built using PAMLX (Yang, 2007) and assuming the HKY85 mutation model (two parameters in the model of DNA evolution) with gamma-distributed rates (approximated by a discrete distribution with eight categories). The Bayesian tree was obtained using BEAST 1.8.3 (Drummond and Rambaut, 2007) and running 50,000,000 iterations, with samples drawn every 10,000 Markov chain Monte Carlo (MCMC) steps. It was visualized using FigTree v.1.4.2. Phylogeny reconstruction was performed considering all the nucleotide substitutions (excluding indels and heteroplasmies) in the coding region (from np 1 to np 14896, relative to mitogenome #1) and the five informative control region mutations 14969, 15159, 15349, 16832A, and 17139A. Thirteen additional control region mutations (see Supplementary information for Figure 1), ten shared exclusively by the previously published mitogenomes #17 and #27 and whose reliability is doubtful, and three (private mutations) seen only in mitogenome #27, were simply superimposed on the magnified MP tree (Figure 1).

The identification of haplogroup diagnostic mutations allowed us to survey 1170 A. albopictus partial COI and ND5 available in the literature for the presence/absence of these mutations. Of these, 284 mtDNAs are from Asia, 349 from the Americas, 32 from Europe, 153 from Africa, and 352 from Oceania (Table 2), thus encompassing populations from both native and non-native areas, and living at different climatic conditions.

The first step in our study was to sequence the entire mtDNA (mitogenome #1) from one mosquito of the Italian (Rimini) laboratory-maintained strain (Table 1; Figure 1), whose nuclear genome was recently sequenced (Dritsou et al., 2015). This mitogenome sequence confirmed that the A. albopictus control region belongs to group 2 of insect control regions (Zhang and Hewitt, 1997), like that of A. aegypti (Dueñas et al., 2006). Three conserved blocks are positioned along the region, which contains two different types (I and II) of tandem repeats (Figure 2). Type I consists of a ∼190 bp repeat unit, whereas type II is made up of a short unit of ∼42 bp. The number of type I and type II tandem repeats that we observed in mitogenome #1 was different from those in sequence NC006817 from Taiwan (mitogenome #27 in Figures 1 and 2) and variable among mtDNAs (Table 1 and data not shown). Moreover, between these two types of repeats, delimited by two conserved blocks, lies an A+T rich region of variable length. The overall length and tandem repeat composition make the PCR amplification and sequencing of the entire A. albopictus control region extremely difficult (see “Material and Methods”). For the reasons outlined above, and difficulties originating from the impossibility of distinguishing heteroplasmy from PCR artifacts (due to replication slippage), we restricted our sequencing of A. albopictus mitogenomes to the coding region (from np 1 to np 14893, NC006817) and nearby control region segments. This approach was employed to obtain the coding region sequences of the additional 24 mitogenomes.

Figure 1 illustrates the Bayesian and MP trees derived from the coding regions of the 27 A. albopictus mitogenomes (25 novel and two previously published). The overall tree structure is virtually identical with the two approaches and is supported by the ML tree (Supplementary Figure 1), indicating a high degree of internal consistency for all major branches. A magnified MP tree is shown in the lower part of Figure 1 in order to illustrate the branch location of the identified mutations.

The mitogenomes (24 distinct haplotypes) cluster into three major branches that we named haplogroups A1, A2, and A3. Haplogroup A1 includes 21 mitogenomes, A2 consists of the five Philippine samples, while A3 encompasses only one mitogenome (#27) from Taiwan. Haplogroups A1 and A2 are rather close to each other and to the A1′2 node. In contrast, the A3 mitogenome differs by 363 mutations from the same node.

Overall, this phylogeny reveals an extensive and previously unreported mitogenome differentiation within A. albopictus. Indeed, when calculated on the standard COI sequence (Hebert et al., 2003) employed for DNA barcoding, the maximum intraspecific divergence was 0.012 (eight mutations in 658 bp), a value in line with those recently reported for A. scutellaris (0.008) and A. aegypti (0.022), but much greater than the value previously reported for A. albopictus (0.002; Sumruayphol et al., 2016).

Haplogroup A1, which encompasses most of the mitogenomes in the phylogeny, is subdivided into two branches that we termed A1a and A1b, with the former further split into A1a1 and A1a2. In our phylogeny these clades and subclades appear to be correlated with different geographic distributions. The branch A1a1 includes the single mitogenome from Japan, the two mitogenomes from the US (Virginia) and many of the mitogenomes from Italy, including one (#1) of the two detected in the Rimini laboratory strain. The mitogenome from Japan (#26) departs from the node A1a1a and its sister clade A1a1a1 contains a sub-branch, A1a1a1a1, of particular interest. It consists of four mitogenomes, two from Italy and the two from the US mentioned above. One of the US mitogenomes (#4) is identical to one from Northern Italy (#5) while the second (#2) is closely related to the mitogenome #3 from Southern Italy.

The sister branch A1a2 is formed by mitogenomes from different regions of Southern Europe (Italy, Albania, and Greece), including the second mitogenome (#14) from the Rimini strain, the previously published Chinese sequence [#17, KR068634 (Zhang et al., 2015)] and the mitogenome (#20) from the Chinese Foshan strain, a laboratory-maintained colony founded in 1981 from mosquitoes from Southeast China. It is worth mentioning that the presence of two distinct haplotypes in different subjects of the Rimini laboratory-maintained strain, one belonging to A1a1a1 and the second to A1a2a1, reveals that at least two females contributed to the genetic formation of the strain. Haplogroup A1b contains all the mitogenomes from Thailand, which cluster in its A1b1 sub-branch, as well as one from Greece and one from Brazil. Finally, haplogroup A2 consists of multiple haplotypes, all from the Philippines.

The phylogenetic analysis not only allowed the identification of A. albopictus haplogroups and sub-haplogroups but also the definition of their distinguishing mutations (Figure 1). These include some diagnostic markers that are located in COI and ND5 partial sequences whose variation has been extensively assessed by published studies and can be retrieved from GenBank (Table 2).

By surveying these sequences for the presence or absence of these mutations, we were able to determine the most likely haplogroup affiliation for most of the 1170 tiger mosquito mtDNAs from populations worldwide. Table 2 reports the haplogroup frequencies for A1a1a1, A1a2, A1b, A2, and A3 obtained from this survey, as well as a category termed “others” that includes mtDNAs that we could not classify and might encompass haplogroups not represented in our phylogeny. Figure 3 provides an overview of the worldwide spatial distribution of these haplogroups.

Phylogenetic analyses revealed that our A. albopictus mitogenomes cluster into the three main haplogroups A1, A2, and A3. Intriguingly, the population screening of the COI mutations (1503, 1578C, 1676C, 1704, 1964) found in mitogenome #27 from Taiwan, the only one belonging to haplogroup A3 in our phylogeny, did not reveal any match in the 1170 tiger mosquito mtDNAs from worldwide populations, not even in the population sample (N = 30) from Taiwan (Table 2). This suggests that A. albopictus mosquitoes with A3 mtDNAs were probably not involved in the recent worldwide spread of the species, raising the possibility that this haplogroup might be rare and/or with a restricted geographical distribution.

In contrast, the survey for COI mutations characterizing the Philippine mitogenomes #21-25 allowed the identification of many other mtDNAs belonging to haplogroup A2 (Table 2), but all in Insular Southeast Asia, suggesting that this haplogroup might be typical and possibly limited to the Philippines, Indonesia, Papua New Guinea, and Northern Australia (Table 2). Therefore, haplogroup A2 appears to have played a role in the spread of A. albopictus from South-East Asia (Figure 3) restricted to the context of Oceania. Indeed A. albopictus is thought to have spread from Indonesia by human-mediated transportation (Beebe et al., 2013). The low frequency of A2 observed on the Papua New Guinean mainland further supports this scenario, whereas the presence of both A1b and A2 mtDNAs along the North Australian border (Table 2) suggests multiple arrivals from distinct geographical sources (Beebe et al., 2013).

Figure 3 shows that, in contrast to the situation described above, members of the other three Asian haplogroups A1a1, A1a2, and A1b are detected in many adventive populations worldwide. This finding identifies these as the Asian mtDNA lineages mainly associated with the recent global spread.

As for haplogroup A1a1, seen in the Japanese mitogenome #26, it is widely distributed in Italy (Figure 1) and shows high frequencies in Central America and Eastern USA. Haplogroup A1a2 is present with frequencies higher than 50% in Japan, Southern China, Singapore, Hawaii, California and Italy, whereas A1b is fixed or almost fixed in Thailand, Malaysia, the Papua New Guinea mainland, Cameroon, and Brazil, but present at much lower frequencies also in California and Texas (Table 2). Even though these geographical distributions are based on the limited population sampling reported in Table 2, some preliminary conclusions can be drawn.

It appears that the ancestral homeland of haplogroup A1a2 might have been a temperate area, possibly Japan or Northern Asia, rather than the tropical range, in agreement with early allozyme studies (Kambhampati et al., 1991). In contrast haplogroup A1b appears to mainly characterize the tropical belt (Figure 3). This may imply that genetic and physiological traits make populations with A1b most suited to the colonization of tropical areas (Birungi and Munstermann, 2002; Kamgang et al., 2011). In our phylogeny the Brazilian A1b mitogenome (#9) harbors a mutational motif that includes the transition at np 7210 in the ND5 gene, also found in all Brazilian samples retrieved from the literature and already identified as a marker for Brazilian A. albopictus (Birungi and Munstermann, 2002). The same transition was found in two ND5 haplotypes retrieved from the literature, one from Phuntsholing in Southern Bhutan (JQ436953) and one from Chiang Mai in Thailand (JQ436956; Porretta et al., 2012), suggesting a probable route of invasion from Indochina (Goubert et al., 2016). Instead the absence of this transition in the samples from Cameroon (Table 2) suggests that A1b mtDNAs arrived in Cameroon from a different tropical source. Interestingly, the absence of photoperiodic diapause in Brazilian mosquitoes supports their origin in tropical Asia, while the diapause in US populations is in agreement with the scenario of an ancestry in the Asian temperate regions (Mousson et al., 2005; Urbanski et al., 2010).

As for haplogroup A1a1, its distribution could not be fully assessed because of the lack of informative marker mutations. However, the absence of the COI transitions at nps 1536 and 2165 distinguishes the members of its main sub-branch, A1a1a1 (Figure 1), from all other mitogenomes in the phylogeny. The survey of these transitions in published data sets revealed that haplogroup A1a1a1 is the most common in Costa Rica, Panama, Texas, and New Jersey (Figure 3) and widespread in Italy (Table 2). This distribution suggests that A1a1a1 and/or its derivatives A1a1a1a and A1a1a1a1 arose recently in an adventive non-Asian population, probably from an ancestral Japanese source, and reached a high frequency because of genetic drift or founder events. From the earliest adventive non-Asian population(s) they then further spread to other distant regions. Such a possibility is in agreement with some previous observations, in particular with the suggestion that Italian tiger mosquitoes, whose first presence was documented in Northern Italy in 1990 (Sabatini et al., 1990), have a dual origin: a possibly direct Northern American source and a probably indirect (through Albania and Greece) Eastern Asian source, with the former related to the international trade of used tires from the eastern coast of the US (Dalla Pozza et al., 1994; Dritsou et al., 2015).

An origin of haplogroup A1a1a1 and/or its derived sub-branches A1a1a1a and A1a1a1a1 in North America and a subsequent arrival to Italy from the US of multiple haplotypes is supported by (i) the detection of four partial sequences from Texas with the transition at np 1823 (Zhong et al., 2013), which defines sub-branch A1a1a1a1, and (ii) our findings that mitogenomes #4 and #5, the first from Virginia and the second from Italy, are identical, and that mitogenomes #2 and #3, again one from Virginia and one from Italy, are closely related (Figure 1). Finally, the scenario that the ancestral source of A1a1a1 might be a northern temperate area such as Japan is further supported by the presence of the COI transitions 1536 and 2435, which characterize mitogenome #26, in five published mtDNAs from Kyoto (JQ004524; Xu and Fonseca, 2011), and is in agreement with allozyme data that have highlighted genetic links between North American, Italian and Japanese populations (Urbanelli et al., 2000).