# Evolutionary history of the little fire ant Wasmannia auropunctata before global invasion: inferring dispersal patterns, niche requirements and past and present distribution within its native range

L. Chifflet, M. S. Rodriguero, L. A. Calcaterra, Olivier REY, V. A. Confalonieri, P. A. Dinghi, F. B. Baccaro, J. L. P. Souza, P. Follett

► To cite this version:

L. Chifflet, M. S. Rodriguero, L. A. Calcaterra, Olivier REY, V. A. Confalonieri, et al., Evolutionary history of the little fire ant Wasmannia auropunctata before global invasion: inferring dispersal patterns, niche requirements and past and present distribution within its native range. Journal of Evolutionary Biology, 2016, 29 (4), pp.790-809. 10.1111/jeb.12827. hal-02637268

HAL Id: hal-02637268

https://hal.inrae.fr/hal-02637268

Submitted on 12 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
doi: 10.1111/jeb.12827

# Evolutionary history of the little fire ant *Wasmannia auropunctata* before global invasion: inferring dispersal patterns, niche requirements and past and present distribution within its native range

L. CHIFFLET*, M. S. RODRIGUERO*, L. A. CALCATERRA†, O. REY‡§, P. A. DINGHI*, F. B. BACCARO¶, J. L. P. SOUZA**, P. FOLLETT†† & V. A. CONFALONIERI*

*Departamento de Ecología, Genética y Evolución, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires (EGE - FCEN - UBA) and IEGEBA (UBA - CONICET), Ciudad Autónoma de Buenos Aires, Argentina

†Fundación para el Estudio de Especies Invasivas (FUEDEI), Hurlingham, Buenos Aires, Argentina

‡INRA, UMR1062, CBGP, Montpellier, France

§CNRS, USR2936, Station d'Ecologie Expérimentale du CNRS à Moulis, Moulis, France

¶Departamento de Biología, Universidade Federal do Amazonas-UFAM, Coroado I, AM, Brazil

**Programa de Pós-Graduação em Entomologia, Instituto Nacional de Pesquisas da Amazônia-INPA, Manaus, AM, Brazil

††USDA-ARS, U.S. Pacific Basin Agricultural Research Center, Hilo, HI, USA

**Keywords:**

Bayesian phylogeography;
biological invasion;
invasive ant;
mitochondrial DNA;
Pleistocene refugia;
range expansion;
South America;
species distribution model.

**Abstract**

The evolutionary history of invasive species within their native range may involve key processes that allow them to colonize new habitats. Therefore, phylogeographic studies of invasive species within their native ranges are useful to understand invasion biology in an evolutionary context. Here we integrated classical and Bayesian phylogeographic methods using mitochondrial and nuclear DNA markers with a palaeodistribution modelling approach, to infer the phylogeographic history of the invasive ant *Wasmannia auropunctata* across its native distribution in South America. We discuss our results in the context of the recent establishment of this mostly tropical species in the Mediterranean region. Our Bayesian phylogeographic analysis suggests that the common ancestor of the two main clades of *W. auropunctata* occurred in central Brazil during the Pliocene. Clade A would have differentiated northward and clade B southward, followed by a secondary contact beginning about 380 000 years ago in central South America. There were differences in the most suitable habitats among clades when considering three distinct climatic periods, suggesting that genetic differentiation was accompanied by changes in niche requirements, clade A being a tropical lineage and clade B a subtropical and temperate lineage. Only clade B reached more southern latitudes, with a colder climate than that of northern South America. This is concordant with the adaptation of this originally tropical ant species to temperate climates prior to its successful establishment in the Mediterranean region. This study highlights the usefulness of exploring the evolutionary history of invasive species within their native ranges to better understand biological invasions.

*Correspondence*: Lucila Chifflet, Departamento de Ecología, Genética y Evolución, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires (EGE - FCEN - UBA), and IEGEBA (UBA-CONICET), Intendente Güiraldes 2160 (C1428EGA), Pabellón II, Ciudad Universitaria, Ciudad Autónoma de Buenos Aires, Argentina.
Tel.: (54 11) 4576 3300 (Ext. 219); fax: (54 11) 4576 3354;
e-mail: lu_chifflet@hotmail.com

**Introduction**

The accidental introduction of exotic pest species is a matter of increasing concern because it often leads to the decline of native species and homogenization of communities (Vitousek *et al.*, 1997; Clavero & García-Berthou, 2005; Lockwood *et al.*, 2007; Sax & Gaines,
2008). Most studies of invasive species have focused on their introduced range (Elton, 1958; Holway *et al.*, 2002; Williams & Smith, 2007), but the importance of also studying their native range is now recognized (Sakai *et al.*, 2001; Bossdorf *et al.*, 2005). Native areas often constitute biologically heterogeneous habitats occupied by ecologically heterogeneous populations (e.g. Winkler *et al.*, 2008), which are constantly subjected to different environmental pressures. On this basis, invasive species may become successful invaders not only due to the release from natural enemies (i.e. predators, parasites, competitors and pathogens) within the introduced areas (Enemy Release Hypothesis, Keane & Crawley, 2002), but also to key evolutionary processes in the native range (Facon *et al.*, 2006; Bossdorf *et al.*, 2008; Lee & Gelembiuk, 2008).

Among the world's most notorious invaders, ants emerge as a serious threat to ecosystem biodiversity, as well as to agriculture and other human activities (Williams, 1994; Holway et al., 2002; O'Dowd et al., 2003). In particular, southern South America has shown to be a major source of invasive ant species to the rest of the world (McGlynn, 1999; Suarez et al., 2010). This is the case for the imported fire ants *Solenopsis invicta* (Ascunce et al., 2011) and *S. richteri* (Tschinkel, 2006), the Argentine ant *Linepithema humile* (Vogel et al., 2010), the dark rover ant *Brachymyrmex patagonicus* (MacGown et al., 2007), the imported crazy ant *Nylan-deria fulva* (LeBrun et al., 2013), the South American big-headed ant *Pheidole obscurithorax* (Wetterer et al., 2015) and the little fire ant *Wasmannia auropunctata* (Rey et al., 2012). The bulk of information on the biology of invasive ants within their native ranges only comes from two intensively studied species, *S. invicta* and *L. humile* (Suarez et al., 2010). Many research has been carried out on the genetics and population structure of imported fire ants (e.g. Ross & Shoemaker, 1997; Ross et al., 1999; Ahrens et al., 2005; Ascunce et al., 2011) and Argentine ants (e.g. Tsutsui & Case, 2001; Vogel et al., 2009), the role of mutualism in the invasion success of Argentine and fire ants (Tillberg et al., 2007; Wilder et al., 2011) and the competitive role of invasive ants in native communities (e.g. LeB-run et al., 2007; Calcaterra et al., 2008). However, no work dealing with the phylogeographic history of these invasive ants within their native range has previously been reported.

*Wasmannia auropunctata* (Roger, 1863) (Hymenoptera: Formicidae), which is native to the Neotropics, was listed as one of the '100 of the World's Worst Invasive Alien Species' (Lowe *et al.*, 2000). During the last 100 years, *W. auropunctata* has become a pest not only in its native range, for example cocoa plantations in Colombia and Brazil (Posada *et al.*, 1976; Delabie, 1989; de Souza *et al.*, 1998), but also in new areas, such as Caribbean and Pacific islands (including Galapagos, Hawaii, New Caledonia, Solomon islands, Tahiti and

Vanuatu), subtropical Atlantic islands (Bermuda and Bahamas), the United States (Florida) and West Africa (e.g. Gabon and Cameroon) (Wetterer & Porter, 2003). It is considered one of the most harmful exotic ant species in the Pacific islands (Loope & Krushelnycky, 2007). Recently, *W. auropunctata* was also found in the Mediterranean region of Israel, demonstrating its potential to invade nontropical habitats (Vonshak *et al.*, 2010). Strikingly, a tropical species could successfully colonize areas that differ in climatic conditions (e.g. cold winter temperatures and seasonality) with respect to its core habitat.

Contrasting to the other species within the genus, *W. auropunctata* has achieved unusual ecological success within its native range (Longino & Fernandez, 2007). Currently, it is broadly distributed from Mexico to northern Argentina inhabiting diverse environments (Wetterer & Porter, 2003). However, it is not clear whether this species is native to this entire region and how extensive the native Neotropical range of *W. auropunctata* was before being spread by humans (Wetterer & Porter, 2003). Mitochondrial DNA analysis indicates that *W. auropunctata* splits into two main clades (A and B) overlapping in Brazil and Ecuador, with clade A distributed over Central America and northern South America and clade B over southern South America (Mikheyev & Mueller, 2007). Clade A populations invaded most regions, except for West Africa, New Caledonia and Tahiti, which were invaded by clade B populations (Mikheyev & Mueller, 2007; Rey et al., 2012).

*Wasmannia auropunctata* provides an interesting model for studying an invasive species which would have undergone major adaptations within the native range before spreading out. A multidisciplinary study conducted by Rey *et al.* (2012) revealed an invasion scenario in which *W. auropunctata* populations at the southern limit of the native distribution range would have adapted to cold temperatures before dispersal to Israel. The genetic evidence in this study suggests that the source population for the Israeli ants originated from clade B and is putatively located in central-eastern Argentina (surrounding the locality of Zárate in Buenos Aires province). Orivel *et al.* (2009) studied demographic and ecological features of tropical *W. auropunctata* populations in a small portion of the native area in French Guiana. They found the coexistence of ecologically heterogeneous populations of this species: (i) non-dominant populations (characterized by small colonies with low worker densities and a few queens per nest) restricted to naturally disturbed areas (mainly flood-plains) of the primary forest and, conversely, (ii) dominant populations (characterized by large colonies with high worker densities and many queens per nest) in anthropogenic areas (e.g. secondary forests or forest edges along roads) associated with low ant species richness, similar to that observed for invasive populations
(e.g. in New Caledonia). In addition, *W. auropunctata* populations within the native range may have experienced a shift in the reproductive system, from nondominant sexual populations in natural habitats to dominant clonal populations occurring mainly in human-modified habitats (Fournier *et al.*, 2005a; Foucaud *et al.*, 2007, 2009, 2010). Some studies have highlighted the role of human land-use in driving evolutionary changes that could promote invasion, such as adaptation to anthropogenic habitats (Hufbauer *et al.*, 2012; Foucaud *et al.*, 2013). Certainly, it is interesting to explore how these ecological, behavioural and physiological changes may have developed over time during the dispersion and diversification of *W. auropunctata* within its native range.

The aim of this study was to investigate the historical demographic processes that have shaped the distribution of *W. auropunctata* within its native range using information based on mitochondrial and nuclear DNA. By focusing on South American populations, we analysed how the species' geographic range has changed over time and examine if the recent introduction in the Mediterranean region is compatible with its phylogenographic history within the native range. We also performed ecological niche modelling (ENM) to study the environmental variables that may have influenced the distribution of *W. auropunctata* within its native range and to determine central and marginal distribution areas. We tested if the southernmost limit of the species' native distribution (i.e. Argentina) was colonized by natural expansion on an evolutionary time-scale, or more recently, by anthropogenic modification of the habitat. We also investigated the potential of *W. auropunctata* for spreading into new suboptimal areas. Finally, we projected the current ecological niche requirements of the species onto past climatic layers to better understand the possible influence of quaternary climatic fluctuations on its current genetic structure, as well as to detect geographic areas that could have acted as Pleistocene refugia.

# Materials and methods

## Sample collection

Due to the wide native range of *W. auropunctata* and to meet the objectives of this study, our sampling focused on the southern populations of the species, covering part of the distribution of clade A and most of the distribution of clade B. Samples were collected from central and northern Argentina, south-eastern Paraguay, southern Bolivia, Uruguay, Brazil, French Guiana and Dominica (Fig. 1). In each of 94 localities, we sampled 1–5 nests, collecting 10–30 workers and 1–5 queens per nest (Table S2). Specimens were preserved in 96% ethanol for DNA sequencing. Voucher specimens were deposited at the Laboratory of Molecular Phylogenies

and Phylogeography (EGE-FCEN-UBA, Argentina). Geographic coordinates and information on the sampling sites (e.g. type of habitat, altitude) were recorded for each nest. Almost all samples were obtained from modified habitats (cities, forests alongside roads, secondary forests, gardens), although a few workers were taken from primary forests.

## Mitochondrial data

### DNA sequencing

Total genomic DNA was extracted from 154 individuals (one worker per nest) using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Inc., Valencia, CA, USA). A 628-bp region of the Cytochrome Oxidase I (COI) gene was amplified by PCR (polymerase chain reaction), using primers LCO and HCO (Folmer *et al.*, 1994). The PCR mixture was a 50 µL-reaction volume containing 50–100 ng of DNA, 3 mM MgCl₂ (Invitrogen), 1× PCR buffer and 1 unit of Taq polymerase (Invitrogen), 0.1 mM of each dNTP (Promega), 0.5 µL of each primer (10 µM) and water to final volume. Thermal cycling conditions were as follows: denaturation at 95 °C for 3 min, then 37 cycles of denaturation at 94 °C for 1 min, annealing at 48 °C for 1 min and extension at 74 °C for 1 min, followed by a final extension at 74 °C for 10 min. The PCR products were visualized on a 1% agarose gel stained with Gel Red (Biotium, Hayward, CA, USA) and purified using a PCR purification kit (Bioneer). The purified PCR products were sequenced with a 3130-XL Automatic Sequencer (Applied Biosystems) at the Sequencing and Genotyping Service of FCEN-UBA. Sequences were inspected, trimmed if necessary and aligned using ClustalW (Thompson *et al.*, 1994). We also checked that the sequences had no gaps or stop-codons that would alter the reading frame, using MEGA v6 (Tamura *et al.*, 2013). All sequences were deposited at GenBank (see Table S1 for voucher specimens and accession numbers (AN.)). The most closely related species *Wasmannia rochai* was used as an outgroup (GenBank AN. EF459824).

### Phylogenetic analysis and gene tree estimation

Unique haplotypes were determined using DNAsp v5.0 (Librado & Rozas, 2009). Phylogenetic relationships between haplotypes were inferred using both maximum parsimony (MP) and probabilistic (Bayesian Analysis (BA)) criteria. The MP analysis was performed using TNT v1.1 (Goloboff *et al.*, 2008). Nodal supports were assessed by bootstrap resampling (1000 replicates). The BA was implemented in BEAST v1.6.2 (Drummond & Rambaut, 2007; Heled & Drummond, 2010). This software allowed us to co-estimate phylogeny and divergence times of the main clades within *W. auropunctata*. We inferred that the best model of sequence evolution under Akaike's information criterion was TrN+I+G using *jMODELTEST2* (Darriba *et al.*, 2012). To
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Map showing sampling locations and distribution of the haplogroups delimited in the phylogenetic analysis. Each symbol represents a locality or too close locations. Clade B haplogroups were represented by circles, subclade A-I by squares and subclade A-II by diamonds. Pie charts of clade B show proportions of haplogroups in their respective colours of haplotypes network of Fig. 3.</figcaption></figure>
further establish the relationships among haplotypes,
we constructed a haplotype network using tcs v1.21
(Clement et al., 2000) with a 95% connection limit.

Genetic diversity

Number of haplotypes (h), haplotype diversity (H_d) and
nucleotide diversity (π) (Nei, 1987) were determined
for each haplogroup identified in the haplotype net-
work.

To detect possible patterns of geographic structure,
COI genetic distances were mapped on the landscape
using Alleles in Space (AIS) v1.0 (Miller, 2005). This pro-
gram was created to analyse from a geographical point of
view the genetic diversity of species with a continuous
and wide distribution, as is the case for W. auropunctata.
The software can be used without defining the popula-
tions a priori because it only requires sample collection of
specimens in different localities throughout the land-
scape under study. The Delaunay triangulation-based
connectivity network was used to link neighbouring
sampling locations together and to calculate genetic dis-
tances (see Guzmán et al., 2012 for further explana-
tions). The AIS program displays a three-dimensional
graph representing a landscape of these genetic distances,
which are expressed as surface heights. Residual genetic
distances were used to avoid possible effects of isolation
by distance. For the total area covered (2204 km²), we
specified an 89 × 24 grid and 0.5 as a distance-weighting
parameter. To better visualize the AIS height output, we
imported the output file into DIVA-GIS v7.1.7 (Hijmans
et al., 2004) and created a two-dimensional colour hot-
spot map overlaid on the geographic study area, with col-
ours corresponding to heights of genetic distance
between points.

Historical demography

To assess whether *W. auropunctata* haplogroups have
experienced recent population expansions, we calcu-
lated Tajima's *D*<sub>T</sub> (Tajima, 1989) and Fu's *F*<sub>S</sub> (Fu, 1997)
indexes using DNAsp v5.0 (Librado & Rozas, 2009). Sig-
nificant deviations from zero were assessed using
10 000 coalescent simulation replicates based on a
Monte Carlo process (Hudson, 1990) and assuming no
recombination in the data sets.

Possible changes in effective population size were also
assessed using the Bayesian Skyline Plot analysis (BSP,
Drummond *et al.*, 2005) implemented in BEAST v1.6.2
(Drummond & Rambaut, 2007). We used a Markov
Chain Monte Carlo (MCMC) of 100 million generations
with a sampling frequency of 1000, under the assump-
tion of a strict clock because this is an intraspecific
study, and the nucleotide evolution model inferred
with JMODELTEST v2.1.4 (Darriba *et al.*, 2012). Using
TRACER (http://tree.bio.ed.ac.uk/software/tracer/), we
checked that the effective sample size (ESS) for each
parameter estimated was > 200 and discarded the first
10% of the trees as the burn-in period. The same

software was then used to perform the Bayesian skyline
reconstruction.

Bayesian phylogeographic analysis

To infer the geographic origin of *W. auropunctata* and its
spatial and temporal distribution during diversification
in South America, we used a relaxed random walk
(RRW) diffusion model across continuous space (Lemey
*et al.*, 2010) in BEAST v1.7.5 (Drummond *et al.*, 2012).
This model incorporates uncertainty in the tree topol-
ogy and the spatial diffusion process to reconstruct
time-sliced contours from the posterior tree distribution
that represent credibility intervals for locations at any
point in time (Lemey *et al.*, 2010).

Exact latitude and longitude for each haplotype were
recorded at the tips of an unknown phylogeny, and the
unobserved two-dimensional locations along each node
of the posterior distribution of the phylogeny were esti-
mated. We added random noise to identical coordinates
using the jitter option with a parameter of 0.001. We
applied the Cauchy RRW model for the continuous trait
model prior (Lemey *et al.*, 2010) and run six indepen-
dent MCMCs for 200 million generations retaining sam-
ples every 10 000 generations using the BEAST priors
specified above. Convergence was evaluated with TRA-
CER, and the chains with ESS > 200 were retained
after discarding the first 10% (i.e. 20 million) iterations
as burn-in. We combined posterior parameter estimates
across the six runs using LogCombiner v1.7.5 and esti-
mated a maximum clade credibility (MCC) tree using
TreeAnnotator v1.7.5 (Heled & Drummond, 2008). The
MCC tree was projected onto the grid of geographical
coordinates using Phylowood (Landis & Bedford, 2014)
to visualize the phylogeographic reconstruction. To esti-
mate the node ages within the phylogenetic inference,
we used a COI divergence rate of 3.54% My⁻¹ (Papa-
dopoulou *et al.*, 2010) because no fossils are available
for this species.

Microsatellite data

To corroborate if nuclear information supported mtDNA
results, we analysed the genotypes of 151 individuals
(111 workers and 40 queens) from the data set of Rey
*et al.* (2012). Our data set consisted of 68 individuals
from clade A and 83 from clade B, of which 79 samples
were from Argentina, 20 from Brazil, 30 from French
Guiana and 22 from Costa Rica. To avoid any impact of
the nest/colony structure on the results, we only used
two workers and one queen per nest for clonal colonies
and two workers and two queens per nest for sexual
colonies. The genotypes consisted of the 12 microsatellite
loci designed by Fournier *et al.* (2005b). The microsatel-
lite matrix is available in Appendix S2.

We employed the Bayesian clustering method imple-
mented in STRUCTURE v2.3.4 (Pritchard *et al.*, 2000)
to investigate how our data were successively grouped
when the number of predetermined genetic clusters (K) was increased from 1 to 18. If the two main clades (A and B) are not an artefact of analysing only mtDNA, we expect that for K = 2 the clustering agrees with the A and B mitochondrial clades. For higher K values, we also evaluated if the clustering coincided with the haplogroup classification based on mtDNA sequence data. All runs involved 3 000 000 MCMC generations, with a burn-in period of 50 000 using an admixture model and the correlated allele frequencies model described by Falush *et al.* (2003). Three runs were performed for each value of K. The number of K that best fits our data was evaluated by calculating ΔK value (Evanno *et al.*, 2005), using Structure Harvester (Earl & vonHoldt, 2012).

It is worth noting that STRUCTURE analysis can only be applied to sexual populations because it assumes the occurrence of panmixia. This limitation prompted us to use this tool for studying the evolution of *K* value based on past gene flow rather than for defining the number of present populations based on the allele frequencies. Although our data set is mainly composed of clonal populations, for the little fire ant there is evidence that these might have arisen recently from local sexual populations (Foucaud *et al.*, 2007), and hence, clonal genotypes most likely reflect recent past sexual genotypes. Therefore, although we searched for the optimal *K* value, we were actually interested in analysing how individuals were successively clustered from *K* = 2 onwards.

## Ecological niche modelling in the Quaternary: past and current distributions

To estimate *W. europunctata* distribution models across Quaternary climatic fluctuations, we implemented the maximum entropy machine-learning algorithm using MAXENT v3.3.3 (Phillips *et al.*, 2006). This program uses presence-only data and environmental variable layers for the studied area. It generates an estimate of probability of species presence ranging from 0 to 1, where 0 is the lowest and 1 the highest probability. Although it is recognized that species distribution is usually controlled by more complex biological processes such as interactions among species, dispersal history and resource heterogeneity (see Fitzpatrick *et al.*, 2007; for limitations of ENMs), the usefulness of these models is contingent on their appropriate use (Araújo & Peterson, 2012). Indeed, ENMs represent a potentially powerful tool in understanding whole-range geography and ecology of species. Ant distribution is known to be largely limited by climate (Roura-Pascual *et al.*, 2011), and this simple model can provide a useful first approach to explore the distribution of *W. europunctata* within its native range.

Past and current climatic variables were downloaded from WorldClim (Hijmans *et al.*, 2005). Climatic vari-

ables for the Last Glacial Maximum (LGM, ≈21 thousand years (kyr) before present (BP)) were obtained from the Paleoclimate Modeling Intercomparison Project Phase II and for the Last Interglacial (LIG, ≈120 kyr BP) from the study by Otto-Bliesner *et al.* (2006). The Holocene scenario (≈6 kyr BP) was not analysed in this work because previous studies dealing with species distribution and palaeoclimatic variables have not found substantial differences between the Holocene and the present (Werneck *et al.*, 2011, 2012). All bioclimatic variables available were tested for multicolinearity by examining cross-correlations among variables (Pearson correlation coefficient, r) within the geographical space based on the occurrence records. The following 11 environmental variables were selected: Annual Mean Temperature, Mean Diurnal Range (Mean of monthly (max temp - min temp)), Isothermality, Maximum Temperature of Warmest Month, Minimum Temperature of Coldest Month, Temperature Annual Range, Mean Temperature of Wettest Quarter, Mean Temperature of Driest Quarter, Annual Precipitation, Precipitation of Wettest Month, and Precipitation of Driest Quarter. LIG and current layers were at a 30 arc-second (≈1 km²) resolution, while LGM layers were at a 2.5 arc-minute resolution. We trained the model including: a) species occurrence data sets within the native range (i.e. excluding locations where the species was introduced), and b) occurrence data sets only for native clade B populations. Data included records from our GPS (Table S2) and Antweb (http://www.antweb.org/bigMap.do?taxonName=myrmicinae). This analysis requires the clade identity of each sample and clade A could not be analysed separately because of the small, insufficient number of specimens that clearly belonged to this clade. We clipped the bioclim rasters to a polygon containing a reasonably estimated native range of *W. europunctata* and then projected onto this same area to determine the potential distribution areas. We repeated this analysis with LIG- and LMG- based projections.

In MAXENT, we used the default convergence threshold (10-5) and increased maximum iterations to 1000; the program automatically selected regularization values and functions of environmental variables. For each run, we used 75% of the localities to train the model and randomly selected 25% of the localities to test the model. We evaluated model performance using a threshold-independent method based on the area under the curve (AUC) of receiver-operating characteristics curve (ROC) (Swets, 1988). The AUC value ranges from 0.5 (random accuracy) to a maximum value of 1 (perfect discrimination). We transformed the output into a map representing probabilities of occurrence because MAXENT produces a continuous probability.
# Results

## Phylogenetic analysis and mtDNA gene tree

The alignment of 154 COI sequences of *W. auropunctata* resulted in 53 unique haplotypes, with 510 constant, 118 polymorphic and 100 parsimony informative characters. Unweighted MP analysis produced 20 equally most parsimonious trees 373 steps long (Fig. S1). Bayesian analyses recovered a tree (Fig. 2) with a topology similar to that of the MP strict consensus tree. In agreement with Mikheyev & Mueller (2007), two main basal clades (A and B) were retrieved with high bootstrap support (MP) and posterior probabilities (BA). Based on both trees and the haplotype network (Fig. 3), we defined two haplogroups within clade A (i.e. groups of monophyletic haplotypes with high posterior probabilities): A-I and A-II, and eight haplogroups within clade

B: B-I to B-VIII. Despite the overlapping distribution of both clades, we found that clade A is mainly distributed in northern South America down to latitude 26°S in San Francisco do Sul, Brazil (Table S2 and Fig. 1), and clade B in southern South America, down to 34°47'S in Lozano (Buenos Aires province, Argentina) (Table S2 and Fig. 1). Within clade A, haplogroup A-II is more extensively distributed than haplogroup A-I. Indeed, the former includes individuals from Paraná, São Paulo, Rio de Janeiro and Bahia states in Brazil, Cordillera state in southern Paraguay, Santa Cruz de la Sierra state in Bolivia, French Guiana and Dominica Island, whereas the latter is only present in French Guiana and Amazonia (Fig. 1). Within clade B, haplogroup B-I emerges basally with respect to the other haplogroups and is found in north-eastern Brazil (Bahia state); B-III, B-VII and B-VIII are mainly located in southern Brazil (but B-VII is also present in eastern Bolivia), B-II occurs in south-eastern Paraguay and central Argentina, B-IV is mainly found in anthropic habitats in central and north-western Argentina and southern Uruguay (but three individuals from Amazonia bear the haplotype H1, i.e. Viruá, Maracá and Ducke); B-V is mainly present in natural habitats in north-eastern Argentina and south-eastern Paraguay, but also includes individuals from the Atlantic forest in the Brazilian state of Santa Catarina and anthropic habitats in north-western Argentina; and B-VI is mainly present in floodplains in north-eastern Argentina, with some individuals from anthropic habitats in central and north-western Argentina (Fig. 1). Interestingly, the most basal haplotype of the B-VI haplogroup was collected in Amazonia (H 47, Manaquiri, Table S2 and Fig. 4).

All the *W. auropunctata* populations from Argentina (the southernmost limit of the species' native distribution) belong to clade B, indicating that it was the only one that could colonize more southern latitudes within the species' native range (Fig. 1).

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Best-fit Bayesian tree obtained with BEAST. Branch colours represent posterior probabilities with a colour gradient from red to blue (maximum and minimum values, respectively). Bootstrap values > 50 are provided for the major clades and subclades (haplogroups defined in this study). The outgroup node was omitted for more clarity.</figcaption></figure>

## Genetic clustering based on microsatellite data

When the number of clusters was set to two (i.e. *K* = 2), the inferred groups are identical to the mtDNA clades A and B (Fig. 5). Only a few individuals show evidence of belonging to both clades (i.e. site 17 (Una, Brazil); Fig. 5), which could suggest historical gene flow between them. Although both clades split with increasing values of *K*, they remain as separate groups. For *K* = 3, clade B populations are grouped in one cluster, whereas clade A is separated into 2 clusters: (i) Costa Rica samples and (ii) Brazil + French Guiana samples. For *K* = 4, both clades are internally subdivided into two groups: clade B is separated into (i) northern Argentina + Itabuna samples, and (ii) central Argentina samples; whereas clade A shows the same subdivision as for *K* = 3. For *K* = 5, clade B shows the same subdivision.
**Fig. 3** Haplotype network constructed with *tcs* 1.21. Each grey circle represents a haplotype (identified by the same number as in Table S1), and circle size is proportional to haplotype frequency. Small white circles represent unique mutation events. Box colours correspond to the different haplogroups.

<figure><img src="image_4.png" /></figure>

**Fig. 4** Result of AIS projected onto the map using DIVA-GIS. Areas in blue and in yellow show the highest and lowest genetic diversity areas, respectively. On the left: distribution of genetic diversity across the sampled area considering both clades A and B. On the right: distribution of genetic diversity considering only clade B.

<figure><img src="image_5.png" /></figure>

vision as for *K* = 4, whereas clade A is internally separated into three groups: (i) Brazil, (ii) French Guiana and (iii) Costa Rica.

The probability of the LnP(D) data increases progressively with increases in *K* in an asymptotic way (Fig. S3). However, *K* = 11 is the most optimal value inferred by Δ*K* (Fig. S4), for which there are no individuals with a mixing of gene pools between clades A and B. The clustering within clade B comprises 5 groups: 2 in northern Argentina, 2 in central Argentina (Buenos Aires) and one in Itabuna (OCR), whereas clade A shows a complex admixture pattern, with many

individuals not assigned to a single population. Nonetheless, we could identify 6 clusters: 2 in Brazil, 2 in French Guiana and 2 in Costa Rica (Fig. 5, *K* = 11).

## Genetic diversity

For the whole data set, haplotype diversity (*H*<sub>d</sub>) was 0.913 and nucleotide diversity (*π*) was 0.03288. Summing nucleotide diversity across clades (Table 1) revealed that clade A has twice the nucleotide diversity of clade B and a slightly higher haplotype diversity, despite the smaller number of samples. Within clade A,
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>
**Fig. 5** Individual grouping inferred from STRUCTURE analysis. Each individual genome is represented by a vertical bar composed of colour sections representing the proportion of membership of each sample to a given population. Only *K* values from 2 to 5 and the most optimal *K* value are shown. At the bottom of the graph, coloured circles and black diamonds and squares next to locality names represent the haplogroup to which each sample belongs according to the mtDNA analysis.

haplogroup A-II has higher nucleotide diversity than A-I despite of its lower haplotype diversity. Within clade B haplogroup B-III, which is endemic to the Atlantic forest, shows the highest genetic diversity. Haplogroup B-VI also exhibits high genetic diversity, probably due to the genetic distance between haplotype 47 from Amazonia and the rest of the haplotypes from Argentina (H4, 5, 12 and 46) (Table S2; Figs 2 and 3). The remaining haplogroups have low genetic diversity values.

The genetic landscape interpolation analysis showed that the genetic diversity of *W. europunctata* within the area under study is higher in the Atlantic forest and northern South America (Fig. 4 left, blue area), and lower in central South America and northern Argentina (Fig. 4 left, yellow area), although there is a small area of high genetic diversity in north-western Argentina. Within clade B, the genetic diversity is higher in northern Argentina, Paraguay, and the Atlantic forest, and lower in central South America and central Argentina (Fig. 4, right).

## Historical demography

Two of the 10 previously defined haplogroups showed signs of population expansion: haplogroups B-IV and B-V (Table 1). Haplogroup B-IV, which is composed of one central and frequent haplotype (H1) with many

<table><caption>**Table 1** Diversity indexes: number of samples (N), number of haplotypes (h), haplotype diversity (H_d), nucleotide diversity (π) calculated for each haplogroup identified in the Bayesian analysis and haplotype network. Fu's F_S and Tajima's D_T indexes were calculated for the haplogroups that presented more than 4 haplotypes.</caption><thead><tr><th>Haplogroup</th><th>N</th><th>h</th><th>H_d</th><th>π</th><th>Fu's F_S</th><th>D_T</th></tr></thead><tbody><tr><td>Clade A</td><td>31</td><td>17</td><td>0.923</td><td>0.02232</td><td>0.127</td><td>-0.19001</td></tr><tr><td>A-I</td><td>5</td><td>5</td><td>1</td><td>0.00796</td><td>-1.345</td><td>-0.95426</td></tr><tr><td>A-II</td><td>26</td><td>12</td><td>0.889</td><td>0.01869</td><td>1.839</td><td>-0.09986</td></tr><tr><td>Clade B</td><td>123</td><td>36</td><td>0.868</td><td>0.01165</td><td>-8.288</td><td>-1.21548</td></tr><tr><td>B-I</td><td>2</td><td>2</td><td>1</td><td>0.00318</td><td>-</td><td>-</td></tr><tr><td>B-II</td><td>6</td><td>2</td><td>0.333</td><td>0.00053</td><td>-</td><td>-</td></tr><tr><td>B-III</td><td>6</td><td>4</td><td>0.8</td><td>0.00881</td><td>1.373</td><td>-0.17192</td></tr><tr><td>B-IV</td><td>58</td><td>13</td><td>0.498</td><td>0.00134</td><td>-10.945**</td><td>-2.34948**</td></tr><tr><td>B-V</td><td>27</td><td>7</td><td>0.687</td><td>0.00153</td><td>-3.069*</td><td>-1.42429</td></tr><tr><td>B-VI</td><td>18</td><td>5</td><td>0.765</td><td>0.00279</td><td>0.165</td><td>-1.17143</td></tr><tr><td>B-VII</td><td>4</td><td>2</td><td>0.667</td><td>0.00106</td><td>-</td><td>-</td></tr><tr><td>B-VIII</td><td>2</td><td>1</td><td>0</td><td>0</td><td>-</td><td>-</td></tr><tr><td>Total</td><td>154</td><td>53</td><td>0.913</td><td>0.03288</td><td>-2.746</td><td>0.0574</td></tr></tbody></table>

less frequent, recently derived haplotypes (Fig. 3), exhibited the strongest signs of population expansion. Haplogroup B-V, which comprises two frequent haplotypes differing by only one mutational step and a few recently derived haplotypes in low frequencies, had a weaker demographic expansion signal. We decided to perform BSP only with individuals of clade B as it has a larger data set and is the only one showing departures from neutrality. The best nucleotide substitution model for this data subset was HKY+I. The BSP for clade B suggested constant population size through time followed by a slight decline and a further rise that seems to continue until present (Fig. 6). Based on the COI divergence rate of 3.54% My⁻¹ (Papadopoulou *et al.*, 2010), the expansion started approximately 28 000 years ago.

## Bayesian phylogeographic analysis

An animation of the inferred distribution changes of *W. europunctata* through time can be viewed in http://mlandis.github.io/phylowood/ by opening the txt matrix (Appendix S1). Bayesian phylogeographic inference suggested that the common ancestor of clades A and B would have lived in central Brazil (Tocantins state) 3.92 million years before present (My BP) (95% highest posterior density 2.67–5.19) (red star in Fig. 7a). A dispersal event from this area would have occurred northward to the Amazon rainforest and southward to the Del Plata river basin and the Atlantic rainforest, giving rise to clades A and B, respectively, between 3.92 and 1.26 My BP (Figs 7a and 8). The common ancestor of clade A was most probably located in the eastern state of Pará in Brazil, at about 1.26 My BP (95% highest posterior density 0.84–1.66). Dispersal and genetic differentiation into multiple adjacent areas possibly took place through the Amazon and San Francisco river basins. The common ancestor of clade B was most probably located further south than the ancestor of clade A, in the north-western state of São Paulo in Brazil, at about 1.22 My BP (95% highest posterior density 0.73–1.69). In Brazil, clade B could have continued its spread and genetic differentiation southward into adjacent areas of eastern Brazil (haplogroup B-I) and southern South America including Argentina, Uruguay, Bolivia and Paraguay (haplogroups B-II – B-VIII) (Figs 7b–e and 8). Note that there exists high uncertainty about the location of the ancestors of each clade represented by the dotted line areas of Fig. 7a. Bayesian analysis suggests a secondary contact between both clades at 0.38–0.16 My BP (Fig. 7d–e) in at least four

<footnote>*P < 0.05; **P < 0.0001.</footnote>
<figure><img src="image_11.png" /></figure>

**Fig. 6** Bayesian Skyline Plot for clade B of *Wasmannia a Europunctata*. The solid line represents the median estimate while the blue area shows the 95% highest posterior density limits. *Y*-axis and *x*-axis show population size and substitutions per site, respectively.

<figure><img src="image_12.png" /></figure>

**Fig. 7** Dispersal patterns of *Wasmannia a Europunctata* within its native range. The common ancestor of all samples and dispersal routes are represented by the red star and lines, respectively. Lines in warmer colours indicate older events and in cooler colours more recent dispersal events. Polygons on the first plot represent 80% highest posterior density regions of the ancestors of clade A and B (i.e. uncertainty about the location of internal nodes in a phylogeny of mtDNA haplotypes). The timing of distribution changes is based on the rate of sequence evolution in insect mtDNA (Papadopoulou *et al.*, 2010), assuming a strict molecular clock.

regions: north-eastern Brazil (Bahia state), southern Brazil (Atlantic rainforest of the Serra do Mar mountain range), Bolivia (Montane forest) and Paraguay (Atlantic semi-deciduous forest). Two back migration events are registered towards northern Brazil (Amazonia) (Fig. 7e), which correspond to the haplotypes H1 (haplogroup B-IV) and H47 (haplogroup B-VI).

## Quaternary palaeomodelling and past and current distributions

According to AUC values, the ability of all models to distinguish presence from random background points was higher than 0.9, indicating good performance (Swets, 1988). The importance of the variables in terms of contribution to the MAXENT model depended on the data set used for training. When using all the native occurrence data set of *W. a Europunctata*, Annual Mean Temperature, Minimum Temperature of Coldest Month and Temperature Annual Range were the most important variables determining species' distribution, with

the former providing the most relevant information. On the other hand, Minimum Temperature of Coldest Month was the most important variable when the data set only included clade B; the other two variables that contributed to the model in decreasing order of importance were Temperature Annual Range and Mean Temperature of Driest Quarter. The response curves of the most important environmental variables in the predictions are shown in Fig. S2a–f. The ecological niche of *W. a Europunctata* (clades A and B) appears to be defined by a Minimum Temperature of the Coldest Month equal or above 23 °C, while that of clade B populations appears to be defined by a Minimum Temperature of the Coldest Month of 10 °C.

The analysis of the potential geographic distribution based on the model trained with all samples indicated that, at present, the areas that were predicted most suitable for the species are in (i) Central America: one spot in Belize and two spots in Nicaragua and Caribbean islands (e.g. Jamaica); and (ii) South America: Bahia state and the Atlantic forest (near São Paulo) in
<figure><img src="image_13.png" /><figcaption>**Fig. 8** Maximum clade credibility tree from Bayesian phylogeographic analysis. Values at nodes indicate node ages (top) and posterior probabilities (bottom). Node bars represent the 95% HPD (highest posterior density) interval and are equivalent to confidence intervals of the node ages. Colour of each haplogroup is congruent with that in the haplotype network of Fig. 3.</figcaption></figure>

Brazil and north-eastern Argentina (Fig. 9, clades A&B). During the LGM, the most probable distribution range of *W. europunctata* was restricted to eastern Mexico and Central America (Costa Rica and Panama) and to eastern Brazil, which corresponds to the Caatinga nuclei, as defined by Prado & Gibbs (1993) or to Pernambuco and Bahia refugia, as defined by Carnaval & Moritz (2008). In addition, a stable area with a lower occurrence probability would have been located in northern Argentina, southern Paraguay and southern Bolivia. Another suitable area was identified in northern South America (northern Amazon block, including parts of Colombia, Venezuela and Guiana) and in southern Brazil. The ENM for the LIG showed a similar distribution pattern to the current period, which is predictable because both of them correspond to interglacial eras.

When the model only included samples belonging to clade B, the distribution pattern became more restricted to southern latitudes (Fig. 9, clade B). At present, the areas with higher probability of clade B occurrence are located in north-western, north-eastern (including the

Mesopotamia) and part of central Argentina (including the mouth of the Paraná river). During the LGM, the species' distribution was restricted to the Atlantic cost of Brazil (the Bahia refugium) and to smaller areas of north-eastern and north-western Argentina (Misiones and the Subandean Piedmont nucleus, as defined by Prado & Gibbs (1993)). Finally, the ENM for the LIG showed a similar distribution pattern to that for the current period.

## Discussion

In this study, we surveyed mitochondrial DNA sequence variation in a 628-bp region of the Cytochrome Oxidase I gene from 154 individuals and nuclear variation at 12 microsatellite loci from 151 individuals of the little fire ant, over a wide region within the species' native range in South America. Our comprehensive data set allowed us to perform a statistical phylogeographic analysis as a contribution to the understanding of the evolutionary history of this highly invasive ant within its native range.
<figure><img src="image_14.png" /><figcaption>**Fig. 9** Potential geographic distribution of *Wasmannia auropunctata* modelled with clade A and B occurrence (top) and clade B occurrence only (bottom), across Quaternary climatic fluctuations: current climatic scenario, during the LGM (21 kyr BP) and the LIG (120 kyr BP). Warmer colours correspond to regions with higher probability of occurrence. White dots show the presence locations used for training, while violet dots show test locations.</figcaption></figure>

**Present distribution and dispersal history of W. auropunctata major clades**

The hypothesis formulated by Mikheyev & Mueller (2007) that there are two genetically distinct clades A and B within *W. auropunctata* was supported by new mitochondrial data generated in this study and greatly strengthened by microsatellite markers. We corroborated that these clades have a disjoint geographic distribution: clade A is almost exclusively tropical, whereas clade B is mostly subtropical and extends into temperate and tropical regions of South America.

Bayesian clustering analysis suggests a strong struc-
turation of our samples, with 11 being the most optimal
*K* value. Such large number of clusters may be attribu-
ted to the great geographic distance between samples,
together with the low dispersal ability of *W. auropunc-
tata*. As previously mentioned, the analysis of the evo-
lution of the *K*-parameter under a Bayesian approach is

useful to explore how individuals are successively
grouped and provides insights into historical gene flow.
For *W. auropunctata*, the hypothesis of separate biologi-
cal units corresponding to two clades was corroborated
because the grouping inferred for *K* = 2 coincided with
the mtDNA-based clade classification. Both major clades
become subdivided with successively increasing *K* val-
ues, but no evidence of recent past gene flow is
observed between them, except for a few individuals
(i.e. population from Una, Brazil).

Bayesian phylogeographic analysis suggests that
*W. auropunctata* originated in central Brazil (with the
highest probability in Tocantins region). The diversifica-
tion into two clades may have occurred during the Plioc-
cene about 3.92 My ago. Two lineages would have
differentiated northward and southward, giving rise to
the ancestors of the current clades A and B, with the
most probable locations being the Amazon and the
Atlantic rainforest, respectively. Several avian species
endemic to the Atlantic forest would have also split from their Amazonian sister species within this same time frame (Cabanne *et al.*, 2008; Miller *et al.*, 2008; Ribas *et al.*, 2009; Patel *et al.*, 2011). These vicariant events can be attributed to the Milankovitch cycles (Bennett, 1990; Haffer, 1993), which led to unstable ecological conditions throughout the Cenozoic, with alternating humid and dry climatic periods associated with continuous vegetation changes (Haffer, 2008). This affected biomes worldwide by isolating populations in different refugia and was an important factor in promoting speciation (Haffer, 1993). The Refugial Speciation Model (Haffer, 1969, 1982, 1993; Müller, 1973; Vanzolini, 1992) was originally applied to tropical forests during the Quaternary and subsequently extended to the Tertiary because it was also affected by Milankovitch cycles (Haffer, 2008). Therefore, the split of several species endemic to the Atlantic and the Amazon forests as well as the separation of *W. auropunctata* into clades A and B could have been driven by these climatic variations. Our study is consistent with the hypothesis that the major diversification events in the Neotropics occurred during the Neogene period (Rull, 2008; Antonelli *et al.*, 2009; Santos *et al.*, 2009). In addition, the establishment of the current Amazon drainage system at approximately 3.0–2.0 My BP (Campbell *et al.*, 2006) could have acted as a geographic barrier contributing to the differentiation between clades A and B of *W. auropunctata*, as indicated by the upper and lower 95% confidence interval for the estimated times of separation between clades. On the other hand, the location of the most common ancestor of the species should be proposed cautiously because the Bayesian phylogeographic inference was not performed with the entire native range of the species; it is possible that results may have changed if more samples from northern and central South America had been included in the analysis. Nonetheless, it is clear that *W. auropunctata* was able to cross the barrier of the Amazon river as indicated by the presence of individuals of the two clades on both sides of the river. The dispersal of *W. auropunctata* across riverine habitats may have been facilitated by the fact that it is well adapted to swamps and floodplains (Foucaud *et al.*, 2009; Orivel *et al.*, 2009).

We detected two sympatric regions for clades A and B; haplogroup A-I coexists with clade B in the Amazonia, whereas haplogroup A-II coexists with clade B in Bolivia, Paraguay and southern Brazil. The Bayesian phylogeographic analysis suggests that the two clades came into secondary contact in these areas at about 380 kyr ago, after they had started geographic expansion and diversification, a process that still continues to the present day. Alternating cycles of cold and warm phases during the Pleistocene (1.8 My BP–10 kyr BP) could have promoted changes in forest composition (Oliveira-Filho & Fontes, 2000) giving rise to these new regions of contact, as hypothesized for

planarians in the Brazilian Atlantic forest (Alvarez-Presas *et al.*, 2014).

Surprisingly, one haplotype of clade B (H1) highly common and widespread in southern South America was also found in northern Brazil (Amazonia). It is difficult to establish whether this is the result of a natural back migration to the north or of a recent human introduction. Likewise, a back migration of clade B to the north seems to be the case for H47, which belongs to haplogroup B-VI and was also found in Amazonia. Further surveys in central Brazil and northern Bolivia are needed to better understand species movements across this region.

## Current most suitable habitats and genetic diversity within *W. auropunctata*

Under the current climatic scenario, ENM including the entire data set (clades A and B) helped us to identify central and marginal distribution areas of the species within the entire native range, on the basis of climatic variables. There are many areas with optimal climatic conditions for the species in Central America, and to a lesser extent, also in Colombia and Venezuela, in north-eastern (Bahia) and south-eastern (São Paulo) Brazil, and north-eastern Argentina (along the Paraná river). As expected, these spots have high genetic diversity, as shown by applying landscape interpolation. Almost the entire Argentinean region is resolved as a marginal area with low genetic diversity. On the other hand, when limiting the data set to only clade B samples, the regions with the most optimal conditions are restricted to southern South America: south-eastern Brazil, north-western (Yungas ecoregion) and north-eastern (Paranaense forest) Argentina and the mouth of the Paraná river. This is in agreement with the southward dispersal of clade B inferred from the Bayesian phylogeographic analysis. Haplogroup B-III (present in south-eastern Brazil) contains the highest genetic diversity within clade B, thus supporting the hypothesis that this clade originated in the Atlantic forest (São Paulo region). The genetic landscape interpolation analysis suggests that the most optimal regions for clade B exhibit a high genetic diversity, except for the mouth of the Paraná river, probably because it was the last colonized area during species expansion. Interestingly, the source population of *W. auropunctata* in Israel would have been located in the vicinity of the mouth of the Paraná river, most likely near the locality of Zárate, Buenos Aires province, Argentina (Rey *et al.*, 2012). This area is dominated by H1, the most common haplotype within clade B, which might be associated with a genetic combination favouring the establishment of *W. auropunctata* in Argentina, the southernmost limit of its native distribution. Here, mean annual temperature is lower than in the rest of the species native range. In this context, it is reasonable to assume that genome evolution at the
southernmost limit of the native distribution of *W. auropunctata* is affected by less benign climate conditions, thus explaining the low genetic diversity in this area.

**Are the southernmost populations of
*W. auropunctata* native or introduced?**

Our study supports that *W. auropunctata* is native to Argentina and was not introduced by human activity based on the following reasons: (i) Bayesian phylogeographic inference of divergence times indicates that the species reached this region between approximately 380–160 kyr BP, and (ii) the high haplotype diversity in Argentina (22 haplotypes) can only be achieved after a long evolutionary process. Clade B is the only one occurring in Argentina, suggesting that this lineage is more adapted to cold and dry climatic conditions than clade A. In fact, the ENM indicates that Minimum Temperature of Coldest Month is the most important environmental variable determining the distribution of the species, and that the probability of the presence of clade B is higher in areas with lower temperatures.

Bayesian Skyline Plot showed that *W. auropunctata* underwent a recent range expansion in Argentina approximately 28 000 years ago. This demographic event could be explained by a global interstadial episode (warmer time period during a glaciation), that took place between 28 000 - 32 000 years ago (Mörner, 1971). Particularly, the values of Tajima's *D*<sub>T</sub> and Fu's *F*<sub>S</sub> suggest demographic expansion for the most distributed haplogroup in Argentina (B-IV). This result, together with genetic landscape interpolation analysis suggesting that genetic diversity is lower in this region than in eastern South America, supports the hypothesis that this marginal area was recently (in evolutionary time) colonized by *W. auropunctata*. Such assumption is in agreement with the star-like pattern for haplogroup B-IV in the haplotype network.

Interestingly, the native range of *W. auropunctata* appears to have expanded southward during the last 50 years (Cuezzo *et al.*, 2015). Indeed, the species was found in regions not predicted by the ENM, with unfavourable climatic conditions, such as the desert of the Monte ecoregion in central-western Argentina (e.g. Anillaco in La Rioja province, Cuezzo *et al.*, 2015). This suggests that the distribution of *W. auropunctata* is probably controlled by nonclimatic factors, meaning that some variables that affect the distribution of the species were not included in the model (e.g. human land-use). *W. auropunctata* could establish in arid regions of western Argentina and, recently, in the Mediterranean region, in the same type of microhabitats (e.g. irrigated gardens and cultivated crops). Hence, human activity may account for the current distribution of *W. auropunctata*, because it is more influenced by temperature than by precipitation, as pointed out by the ENM analysis.

**Possible effect of symbiont infection on mtDNA**

Symbionts are known to play a major role in the genetic structuring and differentiation of the host population. Therefore, results of phylogeographic studies carried out with a single mitochondrial gene from potentially infected populations must be interpreted with caution (Hurst & Jiggins, 2005). In this regard, it is worth mentioning that *Wolbachia* was the only symbiont detected in *W. auropunctata* populations (Rey *et al.*, 2013a). These authors found a strong correlation between infection and reproduction system: most native sexual (nondominant) populations are infected whereas most clonal (dominant) populations are *Wolbachia*-free. Dominant and nondominant populations do not belong to different evolutionary units: trees constructed with either microsatellite or mtDNA data showed no monophyletic clustering, indicating that infected and uninfected populations do not have divergent mtDNA sequences (Foucaud *et al.*, 2007, 2009).

In the present study, most of the sampling was carried out in modified habitats, thus allowing us to assume that most of the populations sampled were uninfected. Indeed, we searched for evidence of *Wolbachia* infection at 25 randomly chosen sites, of which only 6 harboured infected little fire ants (24% of the populations) (data not shown). Furthermore, no haplotype was linked to infection. In particular, samples with the most common haplotype (H1) were not infected with *Wolbachia*. This is in contrast to that expected for an indirect symbiont-mediated selection of a mtDNA variant (see Hurst & Jiggins, 2005). Therefore, these findings together with the congruent clustering obtained with mitochondrial and nuclear markers, avoids any misinterpretation of the results due to symbiont presence in *W. auropunctata* populations.

**Potential historical distribution of *W. auropunctata* based on palaeodistribution modelling**

Our palaeoclimatic modelling agrees with other works supporting the existence of forests in the northern part of the Atlantic forest during the last glaciations (Cabanne *et al.*, 2008) and a recent recolonization in southern Atlantic forest regions (Grazziotin *et al.*, 2006; Cabanne *et al.*, 2007; Carnaval & Moritz, 2008; Carnaval *et al.*, 2009). When both clades were considered, the ENM predicted that the most probable distribution of *W. auropunctata* during the LGM was restricted to eastern Mexico and Central America (Costa Rica and Panama) and to eastern Brazil, corresponding to the Caatinga nucleus (sensu Prado & Gibbs, 1993) and the Atlantic Costal forest. These Brazilian refugia coincide with two other Pleistocene refugia reported in previous works (e.g. Solomon *et al.*, 2008; Carnaval *et al.*, 2009;
Werneck *et al.*, 2011). There seemed to be another stable area for the species in northern Argentina and southern Paraguay and Bolivia (Misiones and Piedmont nuclei (*sensu* Prado & Gibbs, 1993) and Bolivian Interandean forests), in northern South America (northern Amazon block, *sensu* Werneck *et al.* (2011), which includes parts of Colombia, Venezuela and Guiana) and in southern Brazil, which coincides with the São Paulo refugium described by Carnaval *et al.* (2009). The model predicted the presence of a large area of unsuitable habitat across much of the Amazon basin, which was between areas of suitable habitat to the north and south, as is the case for the leafcutter ant *Atta laevigata* (Solomon *et al.*, 2008).

## Evidence for differentiation in niche requirements between clades

For the three climatic periods analysed here, the most suitable regions when considering both clades A and B differed from the ones obtained with clade B only. The former was predicted with high probability in central, northern/north-eastern and southern/south-eastern South America, whereas the latter was predicted only in southern/south-eastern South America. Results suggest that these clades have been subjected and are adapted to different environmental conditions. As it is assumed that niche conservatism between sister taxa breaks down over time (Peterson, 2011), this difference in niche requirement between clades A and B is not surprising taking into account the long time elapsed from the beginning of their differentiation (approximately 3,92 My ago). During the LGM, clade B would have moved northward, with higher, more stable temperatures throughout the year. However, during the interglacial periods, the most suitable areas for this clade seem to have been located at higher southern latitudes supporting the idea that clade B is better adapted to colder conditions than clade A. Some evidence for this statement was provided by experiments showing that workers of clade B populations (from Zárate, Argentina) had significantly higher survival to cold stress compared with workers of clade A populations (from French Guiana) (Rey *et al.*, 2012). Divergent selection or geographic isolation could have led to the diversification of this species. It would be interesting for future studies to contrast differences in behaviour and physiology between these clades, with a focus on genes under positive selection.

## Are clades A and B undergoing a speciation process?

The lack of evidence of gene flow between natural populations of clades A and B, as revealed by microsatellite data, and the large number of mtDNA mutations (21 in total) that separate them, suggest that they could actually constitute different species. Nevertheless,

sexual reproduction between these lineages is possible under laboratory conditions (Rey *et al.*, 2013b), thus indicating that they could be undergoing an incipient speciation process. Further analysis is needed to corroborate this hypothesis. If so, then it is not surprising that an invasive species so far considered typical of tropical climate could establish successfully in a temperate region such as the Mediterranean. The presence of clade B not only in tropical Pacific islands (New Caledonia in 1972 and Tahiti in 1995) and tropical Africa (Gabon in 1914 and Cameroon in 1959), but also in the Mediterranean region (1998) (Foucaud *et al.*, 2010), suggests that it has greater colonization ability than clade A (which invaded only tropical regions), being able to invade a wider spectrum of climate types.

## Taxonomic and biological control implications

The potential origin of *W. auropunctata* in central Brazil is in line with the high diversity of *Wasmannia* species and *W. auropunctata* subspecies in this country (Wetterer & Porter, 2003). It will be interesting to re-evaluate the validity of the different *W. auropunctata* subspecies listed by Kempf (1972), by comparing the correspondence between his classification and the different haplogroups found in our study. Specimens used here are available for further taxonomic and systematic (morphological and genetic) studies to clarify this issue.

The identification of central areas – with the most optimal conditions and the highest genetic diversity – and of the most ancestral areas provides valuable information on the location of natural enemies of this invasive species (Pschorn-Walcher, 1977). To date, the parasitoid *Orasema minutissima* is the only known natural enemy of the little fire ant (Heraty, 1994; Wetterer & Porter, 2003), but its impact on *W. auropunctata* populations has been questioned (Soto *et al.*, 2010). In this regard, our study contributes to the knowledge of areas to search for new candidate biological control agents of *W. auropunctata*.

## Conclusions

This study sheds light on the complex genetic diversity within native populations of the invasive ant *W. auropunctata* and helps to understand the evolutionary processes that shaped its present native distribution. Our results suggest that the two major clades of this species would be evolutionary units which separated approximately 3.92 My ago. ENM under past and current climate scenarios suggests that genetic differentiation between these main clades might have been accompanied by changes in their niche requirements. Clade B, which is the only one extending southward to Argentina, seems to be more adapted to colder temperatures than the exclusively tropical clade A.
It will be worth investigating if the clonal reproduc-
tion system played a key role in the southward expan-
sion of clade B, because clonality has the potential to
fix favourable genetic combinations (e.g. Mergeay *et al.*,
2006). This is likely because almost all the populations
at higher latitudes occurred in anthropogenic habitats,
and many of them were found to be clonal (data not
shown).

The ENM implemented in this work enabled us to
make a comparative analysis with other South Ameri-
can taxa. Our models agree with much previous evi-
dence of different refugia in this region. In addition,
the integrative approach applied here provides new
data at the species level for future testing of hypotheses
that explain the origin of diversity in South America.
Our results are relevant for further comparative studies
concerning phylogeographic patterns for different South
American taxa.

Overall, the present research provides evidence of
how the evolutionary history of *W. a Europunctata* has
influenced its colonization ability and highlights the
possibility that the successful introduction in the
Mediterranean region could have resulted from an
evolutionary process taking place first in its native
range. It also points out that phylogeographic studies
constitute a valuable contribution to better understand
the eco-evolutionary processes underlying invasion
histories.

Acknowledgments

We are grateful to Cristian Tomatis and Noelia Guzmán
for laboratory assistance; to Laura Varone, Priscila
Hanisch and Francisco Sola for providing some of the
samples; and to John Heraty, Edward Lebrun, Dietrich
Gozek and Silvia Pietrokovsky for their helpful com-
ments and suggestions on the manuscript. This research
project was supported by grants from Universidad de
Buenos Aires, ANPCyT and CONICET, Argentina, and
from USDA-ARS-Pacific Basin Agricultural Research
Center. OR is also grateful to the French Ministère de
l'Écologie et du Développement Durable (ECOTROP
program) and the French Agropolis foundation (RTRA-
Montpellier, BIOFIS project 1001-001) for financial
support.

**References**

Ahrens, M.E., Ross, K.G. & Shoemaker, D.D. 2005. Phyloge-
graphic structure of the fire ant *Solenopsis invicta* in its native
South American range: roles of natural barriers and habitat
connectivity. *Evolution* **59**: 1733–1743.

Alvarez-Presas, M., Sánchez-Gracia, A., Carbayo, F., Rozas, J. & Riutort, M. 2014. Insights into the origin and distribution of biodiversity in the Brazilian Atlantic forest hot spot: a statistical phylogeographic study using a low-dispersal organism. *Heredity* **112**: 656–665.

Antonelli, A., Nylanderb, J.A.A., Perssona, C. & Sanmartín, I.
2009. Tracing the impact of the Andean uplift on neotropical
plant evolution. *Proc. Natl. Acad. Sci. USA* **106**: 9749–9754.

Araújo, M.B. & Peterson, A.T. 2012. Uses and misuses of bio-
climatic envelope modelling. *Ecology* **93**: 1527–1539.

Ascunce, M.S., Yang, C.C., Oakey, J., Calcaterra, L.A., Wu,
W.J., Shih, C.J. *et al*. 2011. Global invasion history of the
fire ant *Solenopsis invicta*. *Science* **331**: 1066–1068.

Bennett, K.D. 1990. Milankovitch cycles and their effects on
species in ecological and evolutionary time. *Paleobiology* **16**:
11–21.

Bossdorf, O., Auge, H., Lafuma, L., Rogers, W.E., Siemann, E. & Prati, D. 2005. Phenotypic and genetic differentiation between native and introduced plant populations. *Oecologia* **144**: 1–11.

Bossdorf, O., Lipowsky, A. & Prati, D. 2008. Selection of prea-
dapted populations allowed *Senecio inaequidens* to invade
Central Europe. *Divers. Distrib*. **14**: 676–685.

Cabanne, G.S., Santos, F.R. & Miyaki, C.Y. 2007. Phylogeogra-
phy of *Xiphorhynchus fuscus* (Passeriformes, Dendrocolapti-
dae): vicariance and recent demographic expansion in
southern Atlantic forest. *Biol. J. Linn. Soc.* **91**: 73–84.

Cabanne, G.S., d'Horta, F.M., Sari, E.H.R., Santos, F.R. & Miyaki, C.Y. 2008. Nuclear and mitochondrial phylogeogra-
phy of the Atlantic forest endemic *Xiphorhynchus fuscus* (Aves: Dendrocolaptidae): biogeography and systematics
implications. *Mol. Phylogenet. Evol.* **49**: 760–773.

Calcaterra, L.A., Livore, J.P., Delgado, A. & Briano, J.A. 2008.
Ecological dominance of the red imported fire ant, *Solenopsis
invicta*, in its native range. *Oecologia* **156**: 411–421.

Campbell, K.E. Jr, Frailey, C.D. & Romero-Pittman, L. 2006.
The Pan-Amazonian Ucayali Peneplain, late Neogene sedi-
mentation in Amazonia, and the birth of the modern ama-
zon river system. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **239**:
166–219.

Carnaval, A.C. & Moritz, C.M. 2008. Historical climate change predicts current biodiversity patterns in the Brazilian Atlantic rainforest. *J. Biogeogr.* **35**: 1187–1201.

Carnaval, A.C., Hickerson, M.J., Haddad, C.F.B., Rodrigues,
M.T. & Moritz, C. 2009. Stability predicts genetic diversity in
the Brazilian Atlantic forest hotspot. *Science* **323**: 785–789.

Clavero, M. & García-Berthou, E. 2005. Invasive species are a
leading cause of animal extinctions. *Trends Ecol. Evol.* **20**: 110.

Clement, M., Posada, D. & Crandall, K. 2000. *tcs*: a computer
program to estimate gene genealogies. *Mol. Ecol.* **9**:
1657–1660.

Cuezzo, F., Calcaterra, L.A., Chifflet, L. & Follet, P. 2015. *Was-
mannia* Forel (Hymenoptera: Formicidae: Myrmicinae) in
Argentina: systematics and distribution. *Sociobiology* **62**: 246–
265.

Darriba, D., Taboada, G.L., Doallo, R. & Posada, D. 2012. jMo-
delTest 2: more models, new heuristics and parallel comput-
ing. *Nat. Methods* **9**: 772.

Delabie, J.H.C. 1989. Ocorrência de *Wasmannia a Europunctata*
(Roger 1863) (Hymenoptera, Formicidae, Myrmicinae) em
cacauais na Bahia, Brasil. *Rev. Theob.* **18**: 29–37.

Drummond, A.J. & Rambaut, A. 2007. BEAST: Bayesian evolu-
tionary analysis by sampling trees. *BMC Evol. Biol.* **7**: 214.

Drummond, A.J., Rambaut, A., Shapiro, B. & Pybus, O.G. 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. *Mol. Biol. Evol.* **22**: 1185–1192.
Drummond, A.J., Suchard, M.A., Xie, D. & Rambaut, A. 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. *Mol. Biol. Evol.* **29**: 1969–1973.

Earl, D.A. & vonHoldt, B.M. 2012. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. *Conserv. Genet. Resour.* **4**: 359–361.

Elton, C.S. 1958. *The Ecology of Invasions by Animals and Plants*. Methuen, London.

Evanno, G., Regnaut, S. & Goudet, J. 2005. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. *Mol. Ecol.* **14**: 2611–2620.

Facon, B., Genton, B.J., Shykoff, J., Jarne, P., Estoup, A. & David, P. 2006. A general eco- evolutionary framework for understanding bio- invasions. *Trends Ecol. Evol.* **21**: 130–135.

Falush, D., Stephens, M. & Pritchard, J.K. 2003. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. *Genetics* **164**: 1567–1587.

Fitzpatrick, M.C., Weltzin, J.F., Sanders, N.J. & Dunn, R.R. 2007. The biogeography of prediction error: why does the introduced range of the fire ant over-predict its native range? *Glob. Ecol. Biogeogr.* **16**: 24–33.

Folmer, O., Black, M., Hoeh, W., Lutz, R. & Vrijenhoek, R. 1994. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. *Mol. Mar. Biol. Biotechnol.* **3**: 294–299.

Foucaud, J., Fournier, D., Orivel, J., Delabie, J., Loiseau, A., Le Breton, J.A. *et al.* 2007. Sex and clonality in the little fire ant. *Mol. Biol. Evol.* **24**: 2465–2473.

Foucaud, J., Orivel, J., Fournier, D., Delabie, J.H., Loiseau, A., Le Breton, J. *et al.* 2009. Reproductive system, social organization, human disturbance and ecological dominance in native populations of the little fire ant, *Wasmannia auropunctata*. *Mol. Ecol.* **18**: 5059–5073.

Foucaud, J., Orivel, J., Loiseau, A., Delabie, J.H.C., Jourdan, H., Konghouleux, D. *et al.* 2010. Worldwide invasion by the little fire ant: routes of introduction and eco-evolutionary pathways. *Evol. Appl.* **3**: 363–374.

Foucaud, J., Rey, O., Robert, S., Crespin, L., Orivel, J., Facon, B. *et al.* 2013. Thermotolerance adaptation to human-modified habitats occurs in the native range of the invasive ant *Wasmannia auropunctata* before long-distance dispersal. *Evol. Appl.* **6**: 721–734.

Fournier, D., Estoup, A., Orivel, J., Foucaud, J., Jourdan, H., Le Breton, J. *et al.* 2005a. Clonal reproduction by males and females in the little fire ant. *Nature* **435**: 1230–1235.

Fournier, D., Foucaud, J., Loiseau, A., Cros-Arteil, S., Jourdan, H., Orivel, J. *et al.* 2005b. Characterization and PCR multiplexing of polymorphic microsatellite loci for the invasive ant *Wasmannia auropunctata*. *Mol. Ecol. Notes* **5**: 239–242.

Fu, Y.X. 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. *Genetics* **147**: 915–925.

Goloboff, P.A., Farris, J.S. & Nixon, K.C. 2008. TNT (Tree Analysis Using New Technology) (BETA) ver. 1.1. Published by the authors, Tucumán, Argentina.

Grazziotin, F.G., Monzel, M., Echeverrigaray, S. & Bonatto, S.L. 2006. Phylogeography of the *Bothrops jararaca* complex (Serpentes: Viperidae): past fragmentation and island colonization in the Brazilian Atlantic Forest. *Mol. Ecol.* **15**: 3969–3982.

Guzmán, N.V., Lanteri, A.A. & Confalonieri, V.A. 2012. Colonization ability of two invasive weevils with different reproductive modes. *Evol. Ecol.* **26**: 1371–1390.

Haffer, J. 1969. Speciation in amazonian forest birds. *Science* **165**: 131–137.

Haffer, J. 1982. General aspects of the refuge theory. In: *Biological Diversification in the Tropics* (G.T. Prance, ed), pp. 6–24. Columbia University, New York.

Haffer, J. 1993. Time's cycle and time's arrow in the history of Amazonia. *Biogeographica* **69**: 15–45.

Haffer, J. 2008. Hypotheses to explain the origin of species in Amazonia. *Braz. J. Biol.* **68**: 917–947.

Heled, J. & Drummond, A.J. 2008. Bayesian inference of population size history from multiple loci. *BMC Evol. Biol.* **8**: 289.

Heled, J. & Drummond, A.J. 2010. Bayesian inference of species trees from multilocus data. *Mol. Biol. Evol.* **27**: 570–580.

Heraty, J.M. 1994. Biology and importance of two eucharitid parasites of *Wasmannia* and *Solenopsis*. In: *Exotic Ants: Biology, Impact and Control of Introduced Species* (D. Williams, ed), pp. 104–120. Westview Press, Boulder, CO.

Hijmans, R.J., Guarino, L., Bussink, C., Mathur, P., Cruz, M., Barrentes, I. *et al.* 2004. DIVA-GIS Vsn. 7.1.7. A geographic information system for the analysis of species distribution data. Manual available at: http://www.diva-gis.org.

Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G. & Jarvis, A. 2005. Very high resolution interpolated climate surfaces for global land areas. *Int. J. Climatol.* **25**: 1965–1978.

Holway, D.A., Lach, L., Suarez, A.V., Tsutsui, N.D. & Case, T.J. 2002. The causes and consequences of ant invasions. *Annu. Rev. Ecol. Syst.* **33**: 181–233.

Hudson, R.R. 1990. Gene genealogies and the coalescent process. In: *Oxford Surveys in Evolutionary Biology*, vol. **7** (D. Futuyma & J. Antonovics, eds), pp. 1–44. Oxford University Press, New York.

Hufbauer, R.A., Facon, B., Ravigné, V., Turgeon, J., Foucaud, J., Lee, C. *et al.* 2012. Anthropogenically-Induced Adaptation to Invade (AIAI): contemporary adaptation to human-altered habitats within the native range can promote invasions. *Evol. Appl.* **5**: 89–101.

Hurst, G.D.D. & Jiggins, F.M. 2005. Problems with mitochondrial DNA as a marker in population, phylogeographic, and phylogenetic studies: the effects of inherited symbionts. *Proc. Biol. Sci.* **272**: 1525–1534.

Keane, R.M. & Crawley, M.J. 2002. Exotic plant invasions and the enemy release hypothesis. *Trends Ecol. Evol.* **17**: 164–170.

Kempf, W.W. 1972. Catálogo abreviado das formigas da Região Neotropical. *Stud. Entomol.* **15**: 3–344.

Landis, M.J. & Bedford, T. 2014. Phylowood: interactive web-based animations of biogeographic and phylogeographic histories. *Bioinformatics* **30**: 123–124.

LeBrun, E.G., Tillberg, C.V., Suarez, A.V., Folgarait, P.J., Smith, C.R. & Holway, D.A. 2007. An experimental study of competition between fire ants and Argentine ants in their native range. *Ecology* **88**: 63–75.

LeBrun, E.G., Abbott, J. & Gilbert, L.E. 2013. Imported crazy ant displaces imported fire ant, reduces and homogenizes grassland ant and arthropod assemblages. *Biol. Invasions* **15**: 2429–2442.

Lee, C.E. & Gelembiuk, G.W. 2008. Evolutionary origins of invasive populations. *Evol. Appl.* **1**: 427–448.
Lemey, P., Rambaut, A., Welch, J.J. & Suchard, M.A. 2010. Phylogeography takes a relaxed random walk in continuous space and time. *Mol. Biol. Evol.* **27**: 1877–1885.

Librado, P. & Rozas, J. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. *Bioinformatics* **25**: 1451–1452.

Lockwood, J.L., Hoopes, M.F. & Marchetti, M.P. 2007. *Invasion Ecology*. Blackwell, Oxford.

Longino, J.T. & Fernandez, F. 2007. Taxonomic review of the genus *Wasmannia*. In: *Advances in Ant Systematics (Hymenoptera: Formicidae): Homage to E. O. Wilson – 50 Years of Contributions*, vol. 80 (R.R. Snelling, B.L. Fisher & P.S. Ward, eds), pp. 271–289. Memoires of the American Entomological Institute, Gainesville, FL.

Loope, L.L. & Krushelnycky, P.D. 2007. Current and potential ant impacts in the Pacific Region. *Proc. Hawaii. Entomol. Soc.* **39**: 69–73.

Lowe, S., Browne, M. & Boudjelas, S. 2000. 100 of the world's worst invasive alien species. *Aliens* **12**: 1–12.

MacGown, J.A., Hill, J.G. & Deyrup, M.A. 2007. *Brachymyrmex patagonicus* (Hymenoptera: Formicidae), an emerging pest species in the southeastern United States. *Fla. Entomol.* **90**: 457–464.

McGlynn, T.P. 1999. The worldwide transport of ants: geographic distribution and ecological invasions. *J. Biogeogr.* **26**: 535–548.

Mergaye, J., Verschuren, D. & De Meester, L. 2006. Invasion of an asexual American water flea clone throughout Africa and rapid displacement of native sibling species. *Proc. Biol. Sci.* **273**: 2839–2844.

Mikheyev, A.S. & Mueller, U.G. 2007. Genetic relationships between native and introduced populations of the little fire ant *Wasmannia auropunctata*. *Divers. Distrib.* **13**: 573–579.

Miller, M.P. 2005. Alleles In Space (AIS): computer software for the joint analysis of interindividual spatial and genetic information. *J. Hered.* **96**: 722–724.

Miller, M.J., Bermingham, E., Klicka, J., Escalante, P., do Amaral, F.S., Weir, J.T. et al. 2008. Out of Amazonia again and again: episodic crossing of the Andes promotes diversification in a lowland forest flycatcher. *Proc. Biol. Sci.* **275**: 1133–1142.

Mörner, N.A. 1971. The Position of the Ocean Level during the Interstadial at about 30 000 BP-A Discussion from a Climatic- Glaciologic Point of View. *Can. J. Earth Sci.* **8**: 132–143.

Müller, P. 1973. The dispersal centres of terrestrial vertebrates in the Neotropical realm. *A Study in the Evolution of the Neotropical Biotia and its Native Landscapes*. Junk, The Hague.

Nei, M. 1987. *Molecular Evolutionary Genetics*. Columbia University Press, New York.

O'Dowd, D.J., Green, P.T. & Lake, P.S. 2003. Invasive 'melt-down' on an oceanic island. *Ecol. Lett.* **6**: 812–817.

Oliveira-Filho, A.T. & Fontes, M.A.L. 2000. Patterns of floristic differentiation among Atlantic forests in southeastern Brazil and the influence of climate. *Biotropica* **32**: 793–810.

Orivel, J., Grangier, J., Foucaud, J., Le Breton, J., Andrés, F.X., Jourdan, H. et al. 2009. Ecologically heterogeneous populations of the invasive ant *Wasmannia auropunctata* within its native and introduced ranges. *Ecol. Entomol.* **34**: 504–512.

Otto-Bliesner, B.L., Marshall, S.H., Overpeck, J.T., Miller, G.H. & Hu, A. 2006. Simulating Arctic climate warmth and icefield retreat in the last interglaciation. *Science* **311**: 1751–1753.

Papadopoulou, A., Anastasiou, I. & Vogler, A.P. 2010. Revisiting the insect mitochondrial molecular clock: the mid-Aegean trench calibration. *Mol. Biol. Evol.* **27**: 1659–1672.

Patel, S., Weckstein, J.D., Patane, J.S.L., Bates, J.M. & Aleixo, A. 2011. Temporal and spatial diversification of *Pteroglossus araçaris* (AVES: Ramphastidae) in the Neotropics: constant rate of diversification does not support an increase in radiation during the Pleistocene. *Mol. Phylogenet. Evol.* **58**: 105–115.

Peterson, A.T. 2011. Ecological niche conservatism: a time-structured review of evidence. *J. Biogeogr.* **38**: 817–827.

Phillips, S., Anderson, R.P. & Schapire, R.E. 2006. Maximum entropy modeling of species geographic distributions. *Ecol. Model.* **190**: 231–259.

Posada, L., De Polania, I.Z., De Arevalo, I.S., Saldarriaga, A., Garcia, F. & Cardenas, R. 1976. Lista de insectos dañinos y otras plagas en Colombia. *Instituto Colombiano Agropecuario Boletín Técnico*. **43**: 1–97.

Prado, D.E. & Gibbs, P.E. 1993. Patterns of species distributions in the dry seasonal forests of South America. *Ann. Mo. Bot. Gard.* **80**: 902–927.

Pritchard, J.K., Stephens, P. & Donnelly, P. 2000. Inference of population structure using multilocus genotype data. *Genetics* **155**: 945–959.

Pschorn-Walcher, H. 1977. Biological control of forest insects. *Annu. Rev. Entomol.* **22**: 1–22.

Rey, O., Estoup, A., Vonshak, M., Loiseau, A., Blanchet, S., Calcaterra, L. et al. 2012. Where do adaptive shifts occur during invasion? A multidisciplinary approach to unraveling cold adaptation in a tropical ant species invading the Mediterranean area. *Ecol. Lett.* **15**: 1266–1275.

Rey, O., Estoup, A., Facon, B., Loiseau, A., Aebi, A., Duron, O. et al. 2013a. Distribution of endosymbiotic reproductive manipulators reflects invasion process and not reproductive system polymorphism in the little fire ant *Wasmannia auropunctata*. *PLoS ONE* **8**: e58467.

Rey, O., Facon, B., Foucaud, J., Loiseau, A. & Estoup, A. 2013b. Androgenesis is a maternal trait in the invasive ant *Wasmannia auropunctata*. *Proc. Biol. Sci.* **280**: 20131181.

Ribas, C.C., Miyaki, C.Y. & Cracraft, J. 2009. Phylogenetic relationships, diversification and biogeography in Neotropical Brotogeris parakeets. *J. Biogeogr.* **36**: 1712–1729.

Ross, K.G. & Shoemaker, D.D. 1997. Nuclear and mitochondrial genetic structure in two social forms of the fire ant *Solenopsis invicta*: insights into transitions to an alternate social organization. *Heredity* **78**: 590–602.

Ross, K.G., Shoemaker, D.D., Krieger, M.J., DeHeer, C.J. & Kelner, L. 1999. Assessing genetic structure with multiple classes of molecular markers: a case study involving the introduced fire ant *Solenopsis invicta*. *Mol. Biol. Evol.* **16**: 525–543.

Roura-Pascual, N., Hui, C., Ikeda, T., Ledayd, G., Richardsona, D.M., Carpintero, S. et al. 2011. Relative roles of climatic suitability and anthropogenic influence in determining the pattern of spread in a global invader. *Proc. Natl. Acad. Sci. USA* **108**: 220–225.

Rull, V. 2008. Speciation timing and neotropical biodiversity: the Tertiary – Quaternary debate in the light of molecular phylogenetic evidence. *Mol. Ecol.* **17**: 2722–2729.

Sakai, A.K., Allendorf, F.W., Holt, J.S., Lodge, D.M., Molofsky, J., With, K.A. et al. 2001. The population biology of invasive species. *Annu. Rev. Ecol. Syst.* **32**: 305–332.
Santos, J.C., Coloma, L.A., Summers, K., Caldwell, J.P., Ree, R. & Cannatella, D.C. 2009. Amazonian amphibian diversity is primarily derived from late Miocene Andean lineages. *PLoS Biol.* **7**: e1000056.

Sax, D.F. & Gaines, S.D. 2008. Species invasions and extinction: the future of native biodiversity on islands. *Proc. Natl. Acad. Sci. USA* **105**: 111490–111497.

Solomon, S., Bacci, M., Martins, J., Vinha, G. & Mueller, U. 2008. Paleodistributions and comparative molecular phylo- geography of leafcutter ants (Atta spp.) provide new insight into the origins of Amazonian diversity. *PLoS ONE* **3**: e2738.

Soto, J., Ortiz, J., Murray, E., Morris, B., Varone, L. & Heraty, J.M. 2010. Taxonomy and life history of *Orasema minutissima* (Hymenoptera: Eucharitidae) attacking the little red fire ant in the Caribbean. Entomological Society American Annual Meeting.

de Souza, A.L.B., Delabie, J.H.C. & Fowler, H.G. 1998. *Wasmannia* spp. (Hym., Formicidae) and insect damages to cocoa in Brazilian farms. *J. Appl. Entomol.* **122**: 339–341.

Suarez, A.V., McGlynn, T.P. & Tsuitsui, N.D. 2010. Bioge- graphic and taxonomic patterns of introduced ants. In: *Ant Ecology* (L. Lach, C.L. Parr & K.L. Abbott, eds), pp. 233–244. Oxford University Press, New York.

Swets, J.A. 1988. Measuring the accuracy of diagnostic sys- tems. *Science* **240**: 1285–1293.

Tajima, F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. *Genetics* **123**: 585–595.

Tamura, K., Stecher, G., Peterson, D., Filipski, A. & Kumar, S. 2013. MEGA6: molecular evolutionary genetics analysis version 6.0. *Mol. Biol. Evol.* **30**: 2725–2729.

Thompson, J.D., Higgins, D.G. & Gibson, T.J. 1994. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. *Nucleic Acids Res.* **22**: 4673.

Tillberg, C.V., Holway, D.A., LeBrun, E.G. & Suarez, A.V. 2007. Trophic ecology of invasive Argentine ants in their native and introduced ranges. *Proc. Natl. Acad. Sci. USA* **104**: 20856–20861.

Tschinkel, W.R. 2006. *The Fire Ants*. Harvard University Press, Cambridge.

Tsutsui, N.D. & Case, T.J. 2001. Population genetics and colony structure of the Argentine ant (*Linepithema humile*) in its native and introduced ranges. *Evolution* **55**: 976–985.

Vanzolini, P.E. 1992. Paleoclimas e especiação em animais da América do Sul tropical. *Estud. Av.* **6**: 41–65.

Vitousek, P.M., D'Antonio, C.M., Loope, L.L., Rejmanek, M. & Westbrook, R. 1997. Introduced species: a significant component of human-caused global change. *N. Z. J. Ecol.* **21**: 1–16.

Vogel, V., Pedersen, J.S., d'Ettorre, P., Lehmann, L. & Keller, L. 2009. Dynamics and genetic structure of Argentine ant supercolonies in their native range. *Evolution* **63**: 1627–1639.

Vogel, V., Pedersen, J.S., Giraud, T., Krieger, M.J.B. & Keller, L. 2010. The worldwide expansion of the Argentine ant. *Divers. Distrib.* **16**: 170–186.

Vonshak, M., Dayan, T., Ionescu-Hirsh, A., Freidberg, A. & Hefetz, A. 2010. The little fire ant *Wasmannia auropunctata*: a new invasive species in the Middle East and its impact on the local arthropod fauna. *Biol. Invasions* **12**: 1825–1837.

Werneck, F.P., Costa, G.C., Colli, G.R., Prado, D. & Sites, J.W. Jr 2011. Revisiting the historical distribution of Seasonally

Dry Tropical Forests: new insights based on palaeodistribu- tion modelling and palynological evidence. *Glob. Ecol. Bio- geogr.* **20**: 272–288.

Werneck, F.P., Gamble, T., Colli, G.R., Rodrigues, M.T. & Sites, J.W. Jr 2012. Deep diversification and long-term persistence in the South American 'dry diagonal': integrating continent- wide phylogeography and distribution modeling of geckos. *Evolution* **66**: 3014–3034.

Wetterer, J.K. & Porter, S.D. 2003. The little fire ant *Wasmannia auropunctata*: distribution, impact and control. *Sociobiology* **42**: 1–41.

Wetterer, J.A., MacGown, J.A. & Calcaterra, L.A. 2015. Geo- graphic spread of *Pheidole obscurithorax* (hymenoptera: Formi- cidae). *Trans. Am. Entomol. Soc.* **141**: 222–231.

Wilder, S.M., Holway, D.A., Suarez, A.V. & Eubanks, M.D. 2011. Macronutrient content of plant-based food affects growth of a carnivorous arthropod. *Ecology* **92**: 325–332.

Williams, D.F. (ed.) 1994. *Exotic Ants. Biology, Impact, and Con- trol of Introduced Species*. Westview Press, Boulder, CO.

Williams, S.L. & Smith, J.E. 2007. A global review of the distri- bution, taxonomy, and impacts of introduced seaweeds. *Annu. Rev. Ecol. Evol. Syst.* **38**: 327–359.

Winkler, G., Dodson, J.J. & Lee, C.E. 2008. Heterogeneity within the native range: population genetic analyses of sym- patric invasive and noninvasive clades of the freshwater invading copepod *Eurytemora affinis*. *Mol. Ecol.* **17**: 415–430.

## Supporting information

Additional Supporting Information may be found online in the supporting information tab for this article:

**Figure S1** Maximum parsimony (MP) tree based on *Wasmannia auropunctata* COI haplotypes.

**Figure S2** Response curves of the three most important environmental variables determining the distribution of *W. auropunctata* predicted by the ENM.

**Figure S3** Mean likelihood L(K) and standard deviation per K value obtained by STRUCTURE analysis with three runs per K (2 to 18), using 151 individuals geno- typed for 12 microsatellite loci.

**Figure S4** Evanno *et al.* (2005) ΔK plot for detecting the number of K groups that best fit the microsatellite data of *Wasmannia auropunctata*.

**Table S1** List of sequenced *Wasmannia auropunctata* specimens and GenBank accession numbers.

**Table S2** Sample details: country, locality, coordinates, number of individuals sampled and haplotypes found in each locality.

**Appendix S1** Matrix (in txt. format) to visualize the interactive animation of the Bayesian phylogeographic analysis available at http://mlandis.github.io/phy- lowood/

**Appendix S2** Genotypes of workers (W) and queens (Q) at each of 12 microsatellite loci.

Received 20 April 2015; revised 12 November 2015; accepted 10 January 2016