\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Landscape genetics across the Andes mountains: Environmental variation
drives genetic divergence in the leaf-cutting ant \emph{Atta cephalotes}}



\author[1]{Vanessa Muñoz-Valencia}%
\author[1]{James Montoya-Lerma}%
\author[2]{Perttu Seppä}%
\author[3]{Fernando Diaz}%
\affil[1]{Universidad del Valle}%
\affil[2]{University of Helsinki}%
\affil[3]{Colgate University}%


\vspace{-1em}



  
  \date{May 05, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Distinguishing among the mechanisms underlying the spatial distribution
of genetic variation resulting from the environmental or physical
barriers from those arising due to simple geographic distance is
challenging in complex landscapes. The Andean uplift represents one of
the most heterogeneous habitats where these questions remain unexplored
since multiple mechanisms may interact, confounding their relative
roles. We explore this broad question in the leaf-cutting ant \emph{Atta
cephalotes}, a species that is distributed across the Andes mountains,
using nuclear microsatellite markers and \emph{mtCOI} gene sequences. We
investigate spatial genetic divergence across the western range of the
northern Andes in Colombia by testing the relative role of alternative
scenarios of population divergence, including isolation by geographic
distance (IBD), climatic conditions (IBE), and the physical barriers
presented by the Andes mountains (IBB). Our results reveal substantial
genetic differentiation among \emph{A. cephalotes} populations for both
types of markers, but only nuclear divergence followed a hierarchical
pattern with multiple models of genetic divergence imposed by the
western range. Model selection showed that the IBD, IBE (temperature and
precipitation), and IBB (Andes mountains) models, often proposed as
individual drivers of genetic divergence, interact and explain up to
33\% of the genetic divergence in \emph{A. cephalotes}. The IBE model
remained significant after accounting for IBD, suggesting that
environmental factors play a more prominent role than with IBB. These
factors, in combination with the idiosyncratic dispersal patterns of
ants, appear to determine the hierarchical patterns of gene flow. This
study enriches our understanding of the forces shaping population
divergence in complex habitat landscapes.%
\end{abstract}%



\sloppy


\textbf{Landscape genetics across the Andes mountains: Environmental
variation drives genetic divergence in the leaf-cutting ant \emph{Atta
cephalotes}}

Vanessa Mu\selectlanguage{ngerman}ñoz-Valencia\textsuperscript{1}, James
Montoya-Lerma\textsuperscript{1}, Perttu Seppä\textsuperscript{2} \&
Fernando Diaz\textsuperscript{3}

\textsuperscript{1} Group of Agroecosystem Ecology and Natural Habitats,
Department of Biology, Faculty of Natural Science, Universidad del
Valle, Cali, Colombia.

\textsuperscript{2} Organismal and Evolutionary Biology Research
Programme, Faculty of Biological and Environmental Sciences, University
of Helsinki, Finland.

\textsuperscript{3} Biology Department, Colgate University, Hamilton,
NY, USA.

Corresponding authors:

Vanessa Muñoz-Valencia, vanem28@gmail.com; James Montoya-Lerma,
james.montoya@correounivalle.edu.co; Fernando Diaz, ferdiazfer@gmail.com

\subsection*{ABSTRACT}

{\label{abstract}}

Distinguishing among the mechanisms underlying the spatial distribution
of genetic variation resulting from the environmental or physical
barriers from those arising due to simple geographic distance is
challenging in complex landscapes. The Andean uplift represents one of
the most heterogeneous habitats where these questions remain unexplored
since multiple mechanisms may interact, confounding their relative
roles. We explore this broad question in the leaf-cutting ant \emph{Atta
cephalotes} , a species that is distributed across the Andes mountains,
using nuclear microsatellite markers and \emph{mtCOI} gene sequences. We
investigate spatial genetic divergence across the western range of the
northern Andes in Colombia by testing the relative role of alternative
scenarios of population divergence, including isolation by geographic
distance (IBD), climatic conditions (IBE), and the physical barriers
presented by the Andes mountains (IBB). Our results reveal substantial
genetic differentiation among \emph{A. cephalotes} populations for both
types of markers, but only nuclear divergence followed a hierarchical
pattern with multiple models of genetic divergence imposed by the
western range. Model selection showed that the IBD, IBE (temperature and
precipitation), and IBB (Andes mountains) models, often proposed as
individual drivers of genetic divergence, interact and explain up to
33\% of the genetic divergence in \emph{A. cephalotes} . The IBE model
remained significant after accounting for IBD, suggesting that
environmental factors play a more prominent role than with IBB. These
factors, in combination with the idiosyncratic dispersal patterns of
ants, appear to determine the hierarchical patterns of gene flow. This
study enriches our understanding of the forces shaping population
divergence in complex habitat landscapes.

\subsection*{KEYWORDS}

{\label{keywords}}

Andean uplift, western mountain range, spatial genetic structure,
isolation by distance, isolation by barrier, model selection, isolation
by environment.

\subsection*{INTRODUCTION}

{\label{introduction}}

Population divergence can be determined by the spatial distribution of
individuals, where geographic proximity modulates genetic similarity,
leading to a pattern of isolation by distance -- IBD (Lee \&
Mitchell-Olds, 2011; Shafer \& Wolf, 2013; Slatkin, 1993; Wright, 1943).
However, populations in complex landscapes are often exposed to
environmental variation and physical barriers that can also contribute
to genetic divergence (Manel \& Holderegger, 2013; Manel, Schwartz,
Luikart, \& Taberlet, 2003; Shafer \& Wolf, 2013). Populations could
adapt to their local environment, maximizing the fitness of individuals
under local conditions while decreasing the fitness of immigrants from
alternative environments (Carro, Quintela, Ruiz, \& Barreiro, 2019;
Sobel, 2014; Wang \& Bradburd, 2014; Wang, Glor, \& Losos, 2013). This
adaptive reduction in gene flow can produce a pattern of isolation by
environment -- IBE (Sexton, Hangartner, \& Hoffmann, 2014; Shafer \&
Wolf, 2013; Wang \& Bradburd, 2014). Alternatively, gene flow could be
restricted by allopatric scenarios of genetic differentiation mediated
by physical barriers in the landscape, generating patterns of isolation
by barrier -- IBB (De Queiroz, Torrente-Vilara, Quilodran, da Costa
Doria, \& Montoya-Burgos, 2017; Haffer, 2008; Rull, 2011;
Turchetto-Zolet, Pinheiro, Salgueiro, \& Palma-Silva, 2013). With
increasing landscape complexity, gene flow is likely to be influenced by
a combination of geographical and ecological factors in which these
isolating mechanisms are not mutually exclusive (Crispo, Bentzen,
Reznick, Kinnison, \& Hendry, 2006; Edwards, Keogh, \& Knowles, 2012;
Noguerales, Cordero, \& Ortego, 2016; Wang et al., 2013).

The Andean mountain ranges not only represent one of the most unexplored
environments, but also offer a great complexity of landscapes, promoting
the diversification of a wide range of taxa (Salgado-Roa et al., 2018).
Across these mountains, restricted gene flow has been reported in
multiple organisms, including birds (Cadena, Pedraza, \& Brumfield,
2016), plants (Lagomarsino, Condamine, Antonelli, Mulch, \& Davis, 2016;
Luebert \& Weigend, 2014; Pérez-Escobar et al., 2017), mammals
(Antonelli et al., 2009; Hoorn et al., 2010), insects (Antonelli et al.,
2009; De-Silva et al., 2017; Hoorn et al., 2010), and other arthropods
(Salgado-Roa et al., 2018). However, organisms from contrasting
populations in these habitats are often exposed to a combination of
distance, environmental, and physical barriers to dispersal, challenging
the investigation of the relative roles of different isolating
mechanisms (James, Coltman, Murray, Hamelin, \& Sperling, 2011;
Meirmans, 2015; Noguerales et al., 2016). It remains unclear how the
Andean uplift has modulated patterns of gene flow and the evolution of
several of the most ecologically important groups in the Neotropics,
including social insects. For example, although the entire evolution of
Neotropical ants occurs across the Andes (Mueller et al., 2017), the
interplay between their population structure and environmental variation
relative to the effect of these mountains on isolated populations
remains largely unknown.

The leaf-cutting ant \emph{A. cephalotes} is a major urban and
agricultural pest in the Neotropics, colonizing a wide spectrum of
environments (Della Lucia, Gandra, \& Guedes, 2014; Fernández,
Castro-Huertas, \& Serna, 2015; Hölldobler \& Wilson, 2011). In
Colombia, its distribution overlaps with the maximum complexity of the
Andean uplift, ranging from 0 to 2100 m.a.s.l. (Fernández et al., 2015),
with a vertical thermal gradient of 0.6 °C/100 m (Hermelin, 2015). The
northern section of these mountains in Colombia splits into three main
branches: the western, central, and eastern ranges. Populations
of\emph{A. cephalotes} are separated by these mountains while
simultaneously being exposed to complex combinations of topographical
and environmental variation (Kattan, Franco, Rojas, \& Morales, 2004;
Pérez-Escobar et al., 2017; Salgado-Roa et al., 2018). Such conditions
provide a tremendous climatic spectrum for local adaptation, with an
interplay between population dynamics and species-specific dispersal
patterns (Hakala, Seppä, \& Helanterä, 2019). Evolution under such
environmental heterogeneity could act to shape patterns of gene flow (De
Queiroz et al., 2017; Lee \& Mitchell-Olds, 2011; Noguerales et al.,
2016; Wang et al., 2013), which often produces more complex scenarios
than genetic divergence due to IBD alone (Slatkin, 1993; Wright, 1943).
For example, we recently found that the eastern range of the Andes in
Colombia plays a major role as a geographic barrier to historical gene
flow, restricting the dispersion of \emph{A. cephalotes} from north to
south (Muñoz-Valencia, Vélez-Matínez, Montoya-Lerma, \& Díaz, 2021).
Although this initial study demonstrates the significant influence of
the Andes on population divergence in the leaf-cutting ant at the
phylogeographic scale, the role of local adaptation occurring at more
regional scales across the Andes remains untested.

This study focuses on a finer and more complex geographic distribution
scale of \emph{A. cephalotes} : that of the western range of the Andean
uplift in Colombia. We use a landscape genetic approach to investigate
the role of geographic features and environmental variation in the
definition of patterns of spatial genetic structure in \emph{A.
cephalotes} . Using nuclear (microsatellites) and mitochondrial
(\emph{mtCOI} ) markers, we test the relative roles played by geographic
distance, climate variation, and a major dispersal barrier (the western
range) in modulating patterns of gene flow. As a monogynous
(single-queen) species, \emph{A. cephalotes} presumably undertakes
long-distance nuptial flights that can potentially overcome isolating
barriers (Cherrett, 1968; Helms, 2018; Moser, 1967). Our results
demonstrate that gene flow is limited by a complex interaction of the
three isolating mechanisms (IBD, IBE, and IBB) rather than IBD alone,
while IBE appears to play a stronger role than IBB. Investigating the
spatial genetic structure of a species in an exceptionally heterogeneous
environment helps to elucidate the evolution and diversity of this
ecologically dominant group of ants in the Neotropics.

\subsection*{MATERIALS AND METHODS}

{\label{materials-and-methods}}

\subsubsection*{Sampling}

{\label{sampling}}

Ant sampling was conducted in the Colombian Pacific and Andean regions,
which are separated by the western mountain range of the Colombian Andes
(Figure 1). The Pacific region is classified as a tropical rainforest
with an extremely humid climate and an annual average temperature of 27
°C. The Andean region is further divided into two groups: Andean 1 (800
- 1050 m.a.s.l.) and Andean 2 (1300 - 2200 m.a.s.l.), with highly
variable climatic conditions. The inner valleys in Andean 1 are
climatically classified as tropical savanna and tend to be dry, with an
annual temperature of 25 °C, while the range summits in Andean 2 are
more humid, with a temperate climate, tropical monsoons, and an annual
temperature of 21 °C (Supplementary table ST1) (Chen \& Chen, 2013;
Hernández-Camacho, 1992; Kattan et al., 2004; Köppen, 1884; Peel,
Finlayson, \& McMahon, 2007).

Environmental variation across the three regions was characterized by
differences in temperature, humidity, and precipitation, measured as
five-year averages of the annual temperature (°C), relative humidity
(\%), and precipitation (mm), respectively, for each location (IDEAM,
2019). In addition, a climate classification was represented using four
categories (1 to 4) of different climatic conditions mediated by the
tropical Andes. Variation in topography was estimated by elevation above
sea level. A dummy variable was used to evaluate the Andean uplift as a
major geographic barrier to gene flow. The code 0 was used for
populations from the western side of the western range (Pacific region),
and 1 for populations from the eastern side of these mountains
(Supplementary table ST1).

Worker ants from nine to twenty nests in ten locations (total of 153
nests) were sampled in the period 2017-2018 (Table 1). The distance
between nests in each location was at least 1.5 km, ensuring that the
sampled nests were independent colonies. Three, two, and five locations
were sampled from the Pacific, Andean 1, and Andean 2 regions,
respectively (Table 1).

\subsubsection*{Molecular methods}

{\label{molecular-methods}}

\subsubsection*{\texorpdfstring{\emph{DNA extraction and PCR
amplification of microsatellite
markers}}{DNA extraction and PCR amplification of microsatellite markers}}

{\label{dna-extraction-and-pcr-amplification-of-microsatellite-markers}}

Total DNA extraction was carried out for five workers from each nest
(total 765 workers) using TNES lysis buffer (Tris 50mM, NaCl 0.4M, EDTA
100mM, SDS 0.5\%), pH 7.5, and chloroform:isoamyl alcohol (24:1),
following Wasko et al. (2003), with minor modifications as described by
Muñoz-Valencia et al. (2020).

Thirteen microsatellite loci developed for \emph{A.
cephalotes}(Muñoz-Valencia et al., 2020) were used (Supplementary Table
ST2). PCR reactions were carried out following Muñoz-Valencia et al.
(2020) in a 10 \selectlanguage{greek}µ\selectlanguage{english}L volume containing 10 ng of DNA, 1 X Phusion Flash PCR
Master Mix (Thermo Fisher Scientific), and 2 \selectlanguage{greek}µ\selectlanguage{english}M of each labeled primer.
The thermal profile was: 98 \selectlanguage{ngerman}°C for 1 min, followed by 34 cycles of 98 °C
for 1 s, annealing temperature for 15 s, and 72 °C for 20 s, followed by
a final extension step at 72 °C for 1 min. The fluorescent amplified
fragments were visualized using an automated DNA sequencer ABI 3130
Genetic Analyzer (Applied Biosystems), and allele sizes were estimated
using GeneMapper version 4.0 (Thermo Fisher Scientific).

\subsubsection*{\texorpdfstring{\emph{Sequencing of the mtCOI
gene}}{Sequencing of the mtCOI gene}}

{\label{sequencing-of-the-mtcoi-gene}}

A 368 bp fragment of the mitochondrial cytochrome oxidase subunit I gene
(\emph{mtCOI} ) was sequenced from one sample per nest, obtaining a
total of 146 sequences after discarding failed amplifications and
sequencing (GenBank accession numbers: MW245066 - MW245211). PCR
amplification was carried out using the universal primers Ben and Jerry,
following Kronauer et al. (2004) and Simon et al. (1994). PCR reactions
were performed in a 10 \selectlanguage{greek}µ\selectlanguage{english}l volume, containing 10 ng of DNA, 1X GoTaq\selectlanguage{ngerman}®
Master Mix (PROMEGA), and 2 \selectlanguage{greek}µ\selectlanguage{english}M of each primer. The thermal cycling
profile was 94 \selectlanguage{ngerman}°C for 2 min followed by 30 cycles of 94 °C for 1 min, 58
°C for 1 min, and 72 °C for 1 min, with a final extension step at 72 °C
for 10 min. PCR products were confirmed by electrophoresis on 1 \%
agarose gels. All PCR products were sequenced by Psomagen, Maryland,
USA.

\subsubsection*{Population genetics
analyses}

{\label{population-genetics-analyses}}

\emph{Genetic diversity of DNA microsatellites and the mtCOI gene}

To test for null alleles, large allele dropout, and scoring errors, we
used Microchecker version 2.2.3 (VanOosterhout, Hutchinson, Wills, \&
Shipley, 2004). Putative null alleles were detected in five loci (MAT2,
MAT4, MAT10, MAT15, and MAT28) of five populations (DAG, GUA, DOV, QBY,
and YOB; Table 2), which may have biased the estimates of population
structure (Chapuis \& Estoup, 2007). We accounted for this bias by
comparing analyses after excluding these loci with those of the entire
data set. Our main results and conclusions remained unchanged, however,
regardless of the loci included in the analyses, and therefore only the
results for the entire data set are presented here. Since \emph{A.
cephalotes} nests form family units and nestmates are closely related,
only one individual per nest was used to evaluate the Hardy-Weinberg
equilibrium (HWE) and linkage disequilibrium (LD). Tests were performed
for each locus (HWE) and each locus pair (LD) for each population
(Supplementary Table ST2) using GENEPOP (Rousset, 2008). A Bonferroni
correction was applied with a global significance of \selectlanguage{greek}α\selectlanguage{english}=0.05 for multiple
comparisons.

Genetic diversity was characterized for each locus, region, and
population as allele count (N\textsubscript{a}), allelic richness
(A\textsubscript{R}), expected (H\textsubscript{e}) and observed
(H\textsubscript{o}) heterozygosities, and the inbreeding coefficient
(\emph{F} \textsubscript{IS}), using the software GenAlex v. 6.5
(Peakall and Smouse 2012) and FSTAT v. 2.9.3.2 (Goudet, 1995; 2001)
(Table 2, Supplementary tables ST2, ST3). Tests for significant
differences from zero in the diversity parameters were performed using
the `\emph{aov'}function in R (Chambers, Freeny, \& Heiberger, 1992). We
used the software Bottleneck (Piry et al. 1999) to estimate recent
changes in population sizes by implementing the IAM and SSM mutation
models and a Wilcoxon signed-rank test to evaluate statistical
significances. We also determined whether the allele frequency
distribution differed from the expected L-shape. Where necessary,
significant values were adjusted for multiple comparisons through
Bonferroni correction.

The genetic diversity of \emph{mtCOI} data was estimated through the
number of haplotypes (h), nucleotide diversity (\selectlanguage{greek}π\selectlanguage{english}), and haplotype
diversity (h\textsubscript{d}), using DnaSP 5.19 (Librado \& Rozas,
2009). Tajima's D and Fu's FS neutrality tests were performed to
investigate signatures of recent population expansion (Ramos-Onsins \&
Rozas, 2002) when the null hypothesis of neutrality was rejected due to
significant negative values (\emph{P} \textless{} 0.05 for D,
\emph{P}\textless{} 0.02 for FS; Supplementary table ST3).

\emph{Spatial genetic patterns}

We investigated spatial genetic patterns and restrictions to gene flow
in \emph{A. cephalotes} , using the software ARLEQUIN version 3.5.2.2
(Excoffier \& Lischer, 2010) to perform alternative hierarchical
Analyses of Molecular Variance (AMOVA). First, the populations were
classified into two regions (Pacific and Andean), separated by the
western range of the Colombian Andes, in order to study the effect of
this range as a major geographic barrier to gene flow affecting the
distribution of genetic variation in the data (isolation by barrier,
IBB). Second, the populations were classified into three regions
(Pacific, Andean 1, and Andean 2) defined by the climatic conditions
mediated by the Andes mountains (isolation by environment, IBE;
Supplementary table ST1). In the AMOVA, we estimated the
associated\emph{F} \textsubscript{ST} for microsatellites and the
genetic distance-based \selectlanguage{greek}\emph{Φ} -\selectlanguage{english}statistic for \emph{mtCOI} (Excoffier,
Smouse, \& Quattro, 1992; Meirmans, 2006; Meirmans \& Hedrick, 2011),
both globally and between all pairs of populations. The significance of
the variance components and associated \emph{F\textsubscript{ST}}
and\selectlanguage{greek}\emph{Φ} \selectlanguage{english}indices were calculated using 10,000 permutations.

We further investigated the structure of the \emph{A. cephalotes}nuclear
data using two clustering analyses. First, we used a model-based
Bayesian clustering method implemented in the software STRUCTURE v.
2.3.4 (Pritchard, Stephens, \& Donnelly, 2000), which estimates the
number of genetic clusters (\emph{K} ) independent of spatial sampling.
Analyses were performed using one individual per nest, with and without
admixture, for correlated allele frequencies. A burn-in of 50,000 and
500,000 sampling generations were implemented for \emph{K} ranging from
1 to 12, with 10 iterations for each value of \emph{K} . Evanno's method
(Evanno, Regnaut, \& Goudet, 2005), implemented in STRUCTURE HARVESTER
(Earl \& vonHoldt, 2012), was used to estimate the optimal number of
clusters from the STRUCTURE output (Supplementary figure SF1) and the
results were visualized using CLUMPAK (Kopelman, Mayzel, Jakobsson,
Rosenberg, \& Mayrose, 2015).

Second, discriminant analysis of principal components (DAPC) was
performed using a principal component analysis (PCA) prior to the
discriminant analysis (DA) (Jombart, Devillard, \& Balloux, 2010). The
DA partitions genetic variation, maximizing differences between clusters
while minimizing within-cluster variation. We performed a DAPC analysis
in the R package ADEGENET (Jombart et al., 2010) using one individual
per nest for microsatellite data and a polymorphic nucleotide positions
matrix for \emph{mtCOI} . The function `\emph{dapc} ' was used to
estimate all available principal components (PCs), and to determine the
optimal number of PCs used based on cumulative variance.

\emph{Redundancy analysis: disentangling patterns of isolation by
distance, environment, and barrier}

We implemented a set of redundancy analyses (RDA) to disentangle the
relative contribution of alternative hypotheses for explaining the
genetic differentiation in \emph{A. cephalotes} : IBD, IBE, and IBB. RDA
is a canonical extension of PCA, in which the principal components are
constrained to be linear combinations of a set of predictors (Bradburd,
Ralph, \& Coop, 2013; Legendre \& Fortin, 2010; Meirmans, 2015). The
goal of RDA here was to identify the best ordination model that
describes genetic differentiation (James et al., 2011) in order to
better understand how spatial heterogeneity across the Andes in the
Colombian Pacific and Andean regions affects patterns of gene flow
in\emph{A. cephalotes} .

For the microsatellite data, a conventional RDA analysis was performed
using population allele frequencies as a dependent matrix since the
Euclidean distance estimates involved in the analysis are directly
related to \emph{F} \textsubscript{ST} values (as long as the
frequencies are not scaled) (De Queiroz et al., 2017; Noguerales et al.,
2016). The\emph{mtCOI} data were analyzed through a distance-based RDA,
making direct use of a \selectlanguage{greek}\emph{Φ} \selectlanguage{english}\textsubscript{ST}-based genetic
distance matrix. In both cases, space (geographic distances matrix),
environment (climate variables), and topography (barriers) were
considered explanatory variables. These explanatory variables were
grouped into three classes according to their resulting patterns of
isolation: 1) IBD, with variables representing the geographic distance
between populations (space); 2) IBE, with variables determining
environmental differences between populations; and 3) IBB, based on the
western range as a geographic barrier that splits populations in an
allopatric manner. Since space was the only explanatory variable
initially expressed as a distance matrix, it was transformed into a
vector format (Oksanen et al., 2019) by PCA, using the `\emph{pcnm'}
function in the package VEGAN. Only the best explanatory PCNM components
were retained. The significance of predictors was assessed using
multivariate\emph{F} -statistics with 10000 permutations using the
`\emph{anova.cca} ' function included in the package VEGAN. All
explanatory variables were scaled using the `\emph{scale'} function in
the package VEGAN. The allele frequencies used were not scaled in order
to keep their inter-population relation with the \selectlanguage{greek}\emph{Φ}
\selectlanguage{english}\textsubscript{ST}.

Spatial explanatory variables were tested through an IBD analysis based
on RDA, followed by sequential elimination of environmental variables,
in which only variables with \textbar{} r \textbar{} \textless{} 0.80
were retained (Dormann et al. 2013; Supplementary figure SF2). The
remaining variables were tested through a model selection approach,
comparing all possible combinations or marginal tests to be finally
included in the IBE analysis. Each marginal test was compared to the
null model (intercept; AIC: 1.41). These combinations of environmental
variables were used to identify the best model based on the Akaike
information criteria (AIC). We selected the best model for significant
predictors with the `\emph{ordistep'} function, using the package VEGAN
to prevent overfitting. IBB was tested using the classification of
sampled populations represented by a dummy variable (see sampling). Once
all three models of genetic isolation were defined and tested, both IBE
and IBB were tested while controlling for spatial autocorrelation
through a partial RDA (conditional test). Finally, the relative
contribution of IBD, IBE, and IBB to the explanation of gene flow
patterns in \emph{A. cephalotes} was evaluated from the final complete
model, including all mechanisms of genetic isolation, using a variation
partitioning analysis with the function `\emph{varpart} ' implemented in
the package VEGAN. This analysis allowed us to disentangle total genetic
variance according to the underlying mechanisms of isolation imposed by
space, environment, and barriers, as well as their relative
contributions when considered together.

\subsection*{RESULTS}

{\label{results}}

\subsubsection*{Genetic diversity}

{\label{genetic-diversity}}

All microsatellite loci followed expectations from the HWE or LD tests.
Most loci were highly variable, with the number of alleles ranging from
2 to 20 per locus (Table 2, Supplementary table ST2). The Pacific region
showed the highest allelic richness, largely due to 17 private alleles,
while the Andean 1 region had the lowest allelic richness, with four
private alleles. Most genetic diversity parameters were significantly
higher for the Pacific region, except for the
H\textsubscript{e}estimates (Figure 2A, Supplementary table ST4). An
excess of heterozygotes was detected in the PPY population (IAM model:
T= 3.40,\emph{P} \textless{} 0.01; SSM model: T= 1.92, \emph{P} = 0.03),
while the allele frequencies differed from the L-shaped distribution. A
putative bottleneck was supported by the low number of alleles per locus
detected in this population (range: 1 to 6). Inbreeding coefficients did
not differ significantly from zero, suggesting random mating in all
populations (Table 2).

A 307 bp \emph{mtCOI} gene fragment was successfully sequenced in 146
individuals. Ten nucleotide sites (3.3\%) were polymorphic and defined
10 different haplotypes that differ in 1-3 positions (0.4\% sequence
divergence). The two most common haplotypes, Hap1 (frequency= 0.51) and
Hap4 (frequency = 0.27), were found in almost all sampling locations.
One haplotype was specific to the Pacific (Hap3), one to the Andean 1
(Hap6), and three to the Andean 2 (Hap8, Hap9, Hap10; Figure 2C,
Supplementary Table ST5) regions. Overall, haplotype and nucleotide
diversities were H\textsubscript{d}= 66\% (range: 0.0 -- 76\%) and \selectlanguage{greek}π \selectlanguage{english}=
0.3\% (range: 0.0 -- 0.4\%), respectively. In contrast to the
microsatellite data, the Andean 1 region presented the highest diversity
(Table 2, Figure 2B). However, the differences in mitochondrial
diversity across regions were not significant (Figure 2B, Supplementary
Table ST4).

\subsubsection*{Population structure for nuclear and mitochondrial
data}

{\label{population-structure-for-nuclear-and-mitochondrial-data}}

Two regional models were used in the hierarchical AMOVA, following IBB
and IBE scenarios (Table 3). For the microsatellite data, the results
were similar under IBB and IBE scenarios, showing substantial
differentiation among regions and populations within regions. Moreover,
genetic differentiation was higher among regions when compared with that
among populations, and was slightly higher under the IBE model when
compared to the IBB model (Table 3).

Pairwise \emph{F} \textsubscript{ST} estimates between populations
ranged from 0.01 to 0.28 (mean = 0.10) and were significant following
Bonferroni corrections, except for the BVT-QBD comparison in the Pacific
region (Figure 3A). Pairwise comparisons between regions under the IBE
model were significant, showing the highest differentiation values for
comparisons involving the Andean 1 region (Pacific-Andean 1: 0.04;
Pacific-Andean 2: 0.02; and Andean 1-Andean 2: 0.10), consistent with
the DAPC analysis (Figure 4A). However, regional comparisons under the
IBB model were also significant (\emph{F} \textsubscript{ST} = 0.04).

Hierarchical AMOVA for \emph{mtCOI} data indicated that populations
within regions harbored 34\% of the genetic variance, but no region
effect was detected under the IBB and IBE models (Table 3). The
pairwise\selectlanguage{greek}\emph{Φ} \selectlanguage{english}\textsubscript{ST} values estimated among populations
ranged from 0.00 to 0.86, and 56\% of these comparisons were significant
(Figure 3B). Moreover, DAPC clustering analysis also failed to group
populations by regions (Supplementary figure SF3), supporting the AMOVA
results.

\subsubsection*{Clustering analyses}

{\label{clustering-analyses}}

Evidence for hierarchical population structure in microsatellite data
was also found from the STRUCTURE runs. The \selectlanguage{greek}Δ\selectlanguage{english}\emph{K} \emph{ad
hoc}statistic by Evanno et al. (2005) indicated an optimal \emph{K} = 3
{[}Pr(X / \emph{K} : 2) = -5061.73{]} as the uppermost hierarchical
level (Supplementary figure SF1). These three genetic clusters
correspond to the individuals belonging to the three regions originally
considered under the IBE model, except for the PPY population from the
Andean 2 region. The first group (blue in Figure 4B) clustered most
samples from the Pacific region. The second group (purple in Figure 4B)
clustered most samples from the Andean 1 region, as well as samples from
Andean 2 region (PPY). The third group (orange in Figure 4B) clustered
most samples from Andean 2 region, apart from the PPY population. This
pattern suggests a hierarchical north-south structure in the Andean
region rather than a climate-related structure, as tested through the
IBE-based AMOVA (Figure 4C), and suggests that IBD is influencing
regional differentiation for the Andean regions rather than, or in
addition to, IBE.

To further investigate the role of IBD in the genetic structure
of\emph{A. cephalotes} , we performed a new AMOVA by moving PPY from the
Andean 2 to the Andean 1 region. However, the results were similar to
those produced under the original IBE scenario (results not shown). We
also ran another STRUCTURE analysis, this time only including samples
from the Andean 1 region and the PPY population (southern Andean
cluster). This analysis revealed two genetic groups (K = 2), in which
PPY formed a sole cluster (Figure 4B).

The results for the IBE model were supported by the DAPC analysis, which
clustered defined regions based on environmental variables. The first
two principal components of the DAPC analysis explained 77.1\% of the
variance in allele frequencies (50 PCs retained; Figure 4A). In this
case, the first principal component clustered the Pacific and Andean
regions, while the second principal component reflected the
differentiation between Andean 1 and Andean 2 regions (Figure 4A).

\subsubsection*{RDA analysis -- Alternative scenarios of genetic
differentiation: IBD, IBE, and
IBB}

{\label{rda-analysis-alternative-scenarios-of-genetic-differentiation-ibd-ibe-and-ibb}}

Significant IBD was detected from the RDA analysis for microsatellite
data, in which the geographic distance expressed as the first axis of
the PCNM analysis (PCNM1) contributed to the genetic structure
of\emph{A. cephalotes} , based on model selection. This suggests that
distance is an important determinant of genetic structure in \emph{A.
cephalotes} (Figure 5). These results agree with those obtained from
STRUCTURE, suggesting that genetic differentiation in the Andean regions
can also be explained by an IBD pattern (north to south) and not
exclusively by IBE across regions (Figure 4B). The IBD results followed
different patterns between nuclear and mitochondrial markers since IBD
was not detected from the \emph{mtCOI} data.

Although alternative models used for regional differentiation suggested
the stronger influence of IBE, based on its higher regional genetic
differentiation when compared to IBB or IBD, all of the models were
significant (Table 3). This was also observed from the RDA results,
where two environmental variables (temperature and precipitation) and
the barrier variable (Andes classification) were significantly
associated with genetic divergence (Table 4, Figure 5). Moreover, a
significant contribution of alternative models of genetic
differentiation was only observed from the microsatellite data
(Supplementary table ST6). The optimal model, including all mechanisms
of genetic isolation as tested for IBD, IBE, and IBB, accumulated up to
33\% of explained variation, which remained marginally significant after
accounting for IBD in conditional tests (Figure 5).

Variation partitioning for the full model, including all alternative
models of genetic isolation, showed that IBD and IBE explained the
highest proportion of genetic variation relative to IBB (Figure 5).
Interaction effects between major mechanisms of genetic isolation were
low for the triple interaction (IBD + IBE + IBB = 2\%) and were not
significant for paired interactions (Figure 5). The higher fraction of
unexplained variation (residuals) might reflect the occurrence of
genetic drift within populations, which is not associated with the
explanatory variables.

\subsection*{DISCUSSION}

{\label{discussion}}

The complex landscape across the Andean uplift is an important barrier,
isolating populations and increasing genetic divergence on both sides of
the mountains (Antonelli et al., 2009; Hoorn et al., 2010; Salgado-Roa
et al., 2018). The western range of the Andes in Colombia is one of the
most biodiverse ecosystems in the Neotropics (Kattan et al., 2004;
Salgado-Roa et al., 2018) and is the most complex landscape within the
distribution of the leaf-cutting ant \emph{A. cephalotes} . By
integrating hierarchical population structure with models of isolation
by distance (IBD), environment (IBE), and barrier (IBB), we explored the
role of the western mountain range in the distribution of genetic
variation of \emph{A. cephalotes} . Here, we demonstrated that the
environmental heterogeneity imposed by the Andean uplift has highly
influenced the population structure of \emph{A. cephalotes} .

\subsubsection*{Genetic diversity}

{\label{genetic-diversity-1}}

Our results indicate that the populations of \emph{A. cephalotes} are
highly variable at nuclear markers, with significant genetic
differentiation at both region and population levels. In contrast,
populations are highly differentiated for the \emph{mtCOI} gene, but no
region effect was detected with this marker. Ten \emph{mtCOI} haplotypes
were detected in 146 nests analyzed, showing low nucleotide and moderate
haplotype diversities, with most populations dominated by the same
haplotype (Hap1). As expected considering the common biogeographical
history, similar \emph{mtCOI} patterns have been detected in the
leaf-cutting ant \emph{A. colombica} (Helmkampf, Gadau, \& Feldhaar,
2008), where six \emph{mtCOI} haplotypes and relatively low levels of
nucleotide diversity (0.1\% population-wide) were detected across 20
colonies, with most specimens sharing the same haplotype. However, these
results were obtained from a small sample size in a small area. Despite
the fact that several other genetic studies have been reported in Attini
(\emph{e.g.,} \emph{Acromyrmex} ) (Cantagalli, Mangolin, \&
Ruvolo-Takasusuki, 2013; Diehl, Cavalli-molina, \& Mellender de Araujo,
2002; Diehl, de Ara\selectlanguage{ngerman}újo, \& Cavalli-Molina, 2001; Pinheiro dos Reis,
Fernandes Salomão, de Oliveira Campos, \& Garcia Tavares, 2014), the
nature of the markers and sampling employed prevents us from drawing
meaningful comparisons with these studies regarding diversity estimates.

\subsubsection*{Spatial differentiation at regional and local
levels}

{\label{spatial-differentiation-at-regional-and-local-levels}}

We detected significant genetic differentiation across regions from the
nuclear data, as evidenced from the AMOVA analyses and clustering
methods, suggesting the potential role of the Andean landscape in
restricting gene flow in \emph{A. cephalotes} . We tested alternative
scenarios causing regional differentiation, considering IBE as a model
of divergence, as well as the geographic barrier imposed by the western
range under an IBB model.

Regional differentiation from microsatellite data was higher under the
IBE than the IBB model, suggesting the more prominent role of climate
compared to the geographic barrier in restricting gene flow in \emph{A.
cephalotes} . On the other hand, the STRUCTURE results support a
regional differentiation with mixed effects between IBE and IBD,
evidencing a north-south spatial pattern for Andean populations rather
than a regional climate model. However, this IBD pattern was likely due
to a single population. PPY was initially classified within the Andean 2
region, but clustered mostly (Q \textgreater{} 0.95) with the
geographically closer Andean 1 region. Nevertheless, this pattern occurs
only in the Andean region, suggesting that population differentiation in
this particular case results from dispersal distance rather than
environmental conditions or a geographic barrier imposed by the western
range. Medium to long dispersal distances have been reported for\emph{A.
cephalotes} , with queens flying up to 50 km during nuptial flights
(Cherrett, 1968; Helms, 2018). This distance is at least one order of
magnitude smaller than the inter-population distance used in this study,
but it is also probably an underestimate since it was originally
measured on a small island (Helms, 2018). Moreover, only a small number
of individuals are required in order to homogenize population allele
frequencies, and gene flow can also be conveyed stepwise. Although the
genetic variance at the region level suggested IBE, all of the models of
hierarchical differentiation were significant for the RDA analysis,
suggesting the likelihood that all three leading causes of genetic
isolation are contributing to the observed patterns of gene flow.

Significant nuclear differentiation at the population level suggests
that other features besides proximity, environment, and major dispersal
barriers act to restrict gene flow. The markedly increased variance in
the PPY population is consistent with a recent population expansion
signal. Whether the bottleneck signal is a product of a founder effect
(Barton \& Turelli, 2004; Matute, 2013; Parisod, Trippi, \& Galland,
2005) or a recent change in effective population size cannot be
determined with the current data. However, the proximity and clustering
of PPY to Andean 1 populations rather than the Andean 2 region in the
STRUCTURE results, suggests a recent founder effect or an admixture
between the two Andean regions.

In contrast to the nuclear data, there was no regional differentiation
in the \emph{mtCOI} data, but populations were significantly
differentiated within the regions in both the IBE and IBB models.
Muñoz-Valencia et al. (2021) used \emph{mtCOI} data to study the genetic
structure of \emph{A. cephalotes} from a larger geographic range
covering most of the distribution of the species in South and Central
America. These authors found significant genetic differentiation, both
at the regional scale and among populations within major regions.
Together, these results show that the eastern range of the Andes appears
to be a major dispersal barrier driving regional differentiation, while
the western range forms a relatively homogenous biogeographic area.

\subsubsection*{Disentangling the patterns of IBD, and IBE,
IBB}

{\label{disentangling-the-patterns-of-ibd-and-ibe-ibb}}

Our RDA framework provides strong evidence for restricted gene flow
in\emph{A. cephalotes} across the western mountain range of the
Colombian Andes by all of the mechanisms tested. Furthermore, we found
that the explained genetic variation was maximized when all spatial and
environmental variables associated with these mechanisms were included
in the model. These multifactorial patterns are often expected in
complex and heterogeneous environments (Sexton et al., 2014; Shafer \&
Wolf, 2013), such as in the Colombian Andes (Kattan et al., 2004;
Pérez-Escobar et al., 2017; Salgado-Roa et al., 2018). Regions
classified according to climate variation and geographic barriers follow
IBE but also include IBB, supporting the higher differentiation estimate
(\emph{F} \textsubscript{CT}) obtained under the IBE model and the
significance when both models are included in the RDA (IBE + IBB).
Migration through low elevation passes in the western mountain range
(Hernández-Camacho, 1992; Kattan et al., 2004) would allow gene flow
between regions, explaining why such a major geographic barrier did not
have a more pronounced effect. Spatial and environmental variables
presumably associated with restriction to gene flow often suffer from
spatial autocorrelation, challenging the differentiation of their
underlying causes (Crispo et al., 2006; Edwards et al., 2012; Wang et
al., 2013). Our RDA framework accounted for this issue by demonstrating
strong IBE even after controlling for IBD; while IBB was not
significant, even when ignoring IBD. Temperature and precipitation
therefore appear to be the leading causes of IBE in \emph{A. cephalotes}
.

The significant effect of the environment on the genetic structure
of\emph{A. cephalotes} populations indicates the potential of this
species to adapt to local environmental conditions. This may occur when
processes such as sex-biased dispersal (Edelaar \& Bolnick, 2012) or
selection against migrants (Hendry, 2004; Weber, Bradburd, Stuart,
Stutz, \& Bolnick, 2016; Wright, 1943) decrease the rate of gene flow,
as has been observed in several other species (De Queiroz et al., 2017;
Lee \& Mitchell-Olds, 2011; Wang et al., 2013). Moreover, geographic
distance and barriers in between can restrict gene flow by reducing
dispersal efficiency, which is more associated with local genetic drift
(Clémencet, Viginier, \& Doums, 2005; Cross, Naugle, Carlson, \&
Schwartz, 2016; De Queiroz et al., 2017; Noguerales et al., 2016; Smith
et al., 2018), leading to combined patterns of IBD and IBB. These
complex patterns have been detected, for example, in the Amazonian
common sardine \emph{Triportheus albus} (De Queiroz et al., 2017).
However, disentangling the individual environmental drivers was not
possible due to strong correlations across environmental variables
(Crispo et al., 2006; Edwards et al., 2012; Wang et al., 2013). The
underlying causes of genetic isolation of populations have been rarely
explored in leaf-cutting ants (Branstetter et al., 2017). In \emph{A.
sexdens rubipilosa} , geographic separation of populations alone did not
explain population divergence (Cantagalli et al., 2013), and significant
IBD and biome fragmentation imposed physical barriers to gene flow
in\emph{A. robusta} (Pinheiro dos Reis et al., 2014). All of these
results together strongly suggest that the isolating mechanisms studied
are not mutually exclusive (Crispo et al., 2006; Edwards et al., 2012;
Wang et al., 2013) and are important factors in the evolution of
Neotropical species.

Agricultural activities have significantly contributed to the
colonization processes of leaf-cutting ants (de Carvalho Cabral, 2015).
Conversion of their natural habitat into cultivated lands may favor
population growth and dispersal of leaf-cutting ants (Montoya-Lerma,
Giraldo-Echeverri, Armbrecht, Farji-Brener, \& Calle, 2012; Schowalter
\& Ring, 2017). Queens of \emph{Atta laevigata} have been shown to
prefer nesting in open areas rather than closed forests (Vasconcelos
1997). Our exceptional study population PPY is located in a large city
(Popayán). It is possible that the city infrastructure and constant
human interference, rather than natural dispersal associated with
nuptial flights, have influenced its genetic structure. This pattern is
more compatible with human-mediated dispersal (Zheng, Yang, Zeng, Vargo,
\& Xu, 2018), where environmental changes associated with human
expansion appear to promote population growth in leaf-cutting ants
(Montoya-Lerma et al., 2012; Siqueira et al., 2017). However, PPY is
most likely an exception to a general large-scale migration process
in\emph{A. cephalotes} .

The Andean uplift appears to modulate the population structure
of\emph{A. cephalotes} in a more complex manner than previously thought.
The eastern mountain range of the Andean uplift in Colombia plays a
major role as a geographic barrier to historical gene flow, restricting
the dispersion of \emph{A. cephalotes} from north to south
(Muñoz-Valencia et al., 2021). By exploring a finer scale across the
western mountain range and incorporating neutral genetic markers with
environmental variables, we clearly show that the observed genetic
differentiation in \emph{A. cephalotes} is mainly affected by the
combined multifactorial effect of different isolating mechanisms,
mediated by the landscape complexity. Together, these results suggest
that studies directed at exploring historical gene flow across the Andes
should be interpreted with caution, since the complexity and history of
this landscape can dramatically influence results at different scales.

\subsection*{ACKNOWLEDGEMENTS}

{\label{acknowledgements}}

We thank Kirsi Kähkönen, the ESB (Evolution, Sociality, and Behavior)
research group, and everyone in the Molecular Ecology and Systematics
laboratory (University of Helsinki) for helping to develop the DNA
microsatellite loci. We also thank Sandra Milena Valencia Giraldo,
Andrea López-Peña, and Glever Alexander Vélez-Martínez for helping with
the DNA extraction. This work was funded by the Vice-Rectoría de
Investigaciones, Universidad del Valle, Cali, Colombia (grant number:
CI71067); and COLCIENCIAS National Program of PhD (grant number:
617-2013).

\subsection*{REFERENCES}

{\label{references}}

Antonelli, A., Quijada-Mascareñas, A., Crawford, A. J., Bates, J. M.,
Velazco, P. M., \& Wüster, W. (2009). Molecular studies and
phylogeography of Amazonian tetrapods and their relation to geological
and climatic models. In C. Hoorn \& F. Wesselingh (Eds.),\emph{Amazonia,
landscape and species evolution: A look into the past}(pp. 386--404).
London, UK: Blackwell Publishing Ltd.

Barton, N. H., \& Turelli, M. (2004). Effects of genetic drift on
variance components under a general model of epistasis.\emph{Evolution}
, \emph{58} (10), 2111--2132. doi: 10.1111/j.0014-3820.2004.tb01591.x

Bradburd, G. S., Ralph, P. L., \& Coop, G. M. (2013). Disentangling the
effects of geographic and ecological isolation on genetic
differentiation. \emph{Evolution} , \emph{67} (11), 1--25. doi:
10.1111/evo.12193.Disentangling

Branstetter, M. G., Jesovnik, A., Sosa-Calvo, J., Lloyd, M. W.,
Faircloth, B. C., Brady, S. G., \& Schultz, T. R. (2017). Dry habitats
were crucibles of domestication in the evolution of agriculture in
ants.\emph{Proceedings of the Royal Society B} , \emph{284} , 20170095.

Cadena, C. D., Pedraza, C. A., \& Brumfield, R. T. (2016). Climate,
habitat associations and the potential distributions of Neotropical
birds: Implications for diversification across the Andes. \emph{Revista
de La Academia Colombiana de Ciencias Exactas, Físicas y Naturales}
,\emph{40} (155), 275. doi: 10.18257/raccefyn.280

Cantagalli, L. B., Mangolin, C. A., \& Ruvolo-Takasusuki, M. C. C.
(2013). Population genetics of Atta sexdens rubropilosa (Hymenoptera:
Formicidae). \emph{Acta Biológica Colombiana} , \emph{18} (1), 179--190.

Carro, B., Quintela, M., Ruiz, J. M., \& Barreiro, R. (2019). Wave
exposure as a driver of isolation by environment in the marine gastropod
Nucella lapillus. \emph{Hydrobiologia} , \emph{839} , 51--69. doi:
10.1007/s10750-019-03993-5

Chambers, J. M., Freeny, A., \& Heiberger, R. M. (1992). Analysis of
variance; designed experiments. In J. M. Chambers \& T. J. Hastie
(Eds.), \emph{Statistical models} . doi: 10.1007/978-3-642-50096-1

Chapuis, M. P., \& Estoup, A. (2007). Microsatellite null alleles and
estimation of population differentiation. \emph{Molecular Biology and
Evolution} , \emph{24} (3), 621--631. doi: 10.1093/molbev/msl191

Chen, D., \& Chen, H. W. (2013). Using the Köppen classification to
quantify climate variation and change: An example for
1901-2010.\emph{Environmental Development} , \emph{6} (1), 69--79. doi:
10.1016/j.envdev.2013.03.007

Cherrett, J. (1968). A flight record for queens of Atta cephalotes L.
(Hym. Formicidae). \emph{Entomologist 's Monthly Magazine} , \emph{104}
, 255--256.

Clémencet, J., Viginier, B., \& Doums, C. (2005). Hierarchical analysis
of population genetic structure in the monogynous ant Cataglyphis cursor
using microsatellite and mitochondrial DNA markers. \emph{Molecular
Ecology} , \emph{14} (12), 3735--3744. doi:
10.1111/j.1365-294X.2005.02706.x

Crispo, E., Bentzen, P., Reznick, D. N., Kinnison, M. T., \& Hendry, A.
P. (2006). The relative influence of natural selection and geography on
gene flow in guppies. \emph{Molecular Ecology} , \emph{15} (1), 49--62.
doi: 10.1111/j.1365-294X.2005.02764.x

Cross, T. B., Naugle, D. E., Carlson, J. C., \& Schwartz, M. K. (2016).
Hierarchical population structure in greater sage-grouse provides
insight into management boundary delineation. \emph{Conservation
Genetics} , \emph{17} (6), 1417--1433. doi: 10.1007/s10592-016-0872-z

De-Silva, D. L., Mota, L. L., Chazot, N., Mallarino, R., Silva-Brandão,
K. L., Gómez-Piñerez, L. M., \ldots{} Elias, M. (2017). North Andean
origin and diversification of the largest ithomiine butterfly
genus.\emph{Scientific Reports} , \emph{7} , 1--17. doi:
10.1038/srep45966

de Carvalho Cabral, D. (2015). Into the bowels of tropical earth:
Leaf-cutting ants and the colonial making of agrarian
Brazil.\emph{Journal of Historical Geography} , \emph{50} , 92--105.
doi: 10.1016/j.jhg.2015.06.014

De Queiroz, L. J., Torrente-Vilara, G., Quilodran, C., da Costa Doria,
C. R., \& Montoya-Burgos, J. I. (2017). Multifactorial genetic
divergence processes drive the onset of speciation in an Amazonian
fish.\emph{PLoS ONE} , \emph{12} (12), 1--27. doi:
10.1371/journal.pone.0189349

Della Lucia, T. M., Gandra, L. C., \& Guedes, R. N. (2014). Managing
leaf-cutting ants: Peculiarities, trends and challenges. \emph{Pest
Management Science} , \emph{70} (1), 14--23. doi: 10.1002/ps.3660

Diehl, E., Cavalli-molina, S., \& Mellender de Araujo, A. (2002).
Isoenzyme variation in the leaf-cutting ants Acromyrmex heyeri and
Acromyrmex striatus (Hymenoptera, Formicidae). \emph{Genetics and
Molecular Biology} , \emph{25} (2), 173--178.

Diehl, E., de Araújo, A. M., \& Cavalli-Molina, S. (2001). Genetic
variability and social structure of Colonies in Acromyrmex heyeri and A.
striatus (Hymenoptera: Formicidae). \emph{Brazilian Journal of Biology =
Revista Brasleira de Biologia} , \emph{61} (4), 667--678. doi:
10.1590/S1519-69842001000400017

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré,
G., \ldots{} Lautenbach, S. (2013). Collinearity: A review of methods to
deal with it and a simulation study evaluating their
performance.\emph{Ecography} , \emph{36} (1), 027--046. doi:
10.1111/j.1600-0587.2012.07348.x

Earl, D. A., \& vonHoldt, B. M. (2012). STRUCTURE HARVESTER: a website
and program for visualizing STRUCTURE output and implementing the Evanno
method. \emph{Conservation Genetics Resources} , \emph{4} (2), 359--361.
doi: 10.1007/s12686-011-9548-7.

Edelaar, P., \& Bolnick, D. I. (2012). Non-random gene flow: An
underappreciated force in evolution and ecology. \emph{Trends in Ecology
and Evolution} , \emph{27} (12), 659--665. doi:
10.1016/j.tree.2012.07.009

Edwards, D. L., Keogh, J. S., \& Knowles, L. L. (2012). Effects of
vicariant barriers, habitat stability, population isolation and
environmental features on species divergence in the south-western
Australian coastal reptile community. \emph{Molecular Ecology}
,\emph{21} (15), 3809--3822. doi: 10.1111/j.1365-294X.2012.05637.x

Evanno, G., Regnaut, S., \& Goudet, J. (2005). Detecting the number of
clusters of individuals using the software STRUCTURE: A simulation
study. \emph{Molecular Ecology} , \emph{14} (8), 2611--2620. doi:
10.1111/j.1365-294X.2005.02553.x

Excoffier, L., \& Lischer, H. E. L. (2010). Arlequin suite ver 3.5: A
new series of programs to perform population genetics analyses under
Linux and Windows. \emph{Molecular Ecology Resources} , \emph{10} (3),
564--567. doi: 10.1111/j.1755-0998.2010.02847.x

Excoffier, L., Smouse, P. E., \& Quattro, J. M. (1992). Analysis of
molecular variance inferred from metric distances among DNA haplotypes:
Application to human mitochondiral DNA restriction data.\emph{Genetics}
, \emph{131} , 479--491. doi: 10.1007/s00424-009-0730-7

Fernández, F., Castro-Huertas, V., \& Serna, F. (2015). Hormigas
cortadoras de hojas de Colombia: Acromyrmex \& Atta (Hymenoptera:
Formicidae). In \emph{Fauna de Colombia} (Monografía). Bogotá: Instituto
de Ciencias Naturales, Universidad Nacional de Colombia.

Goudet, J. (1995). FSTAT (Version 1.2): A Computer program to Calculate
F-Statistics. \emph{Journal of Heredity} , pp. 485--486.

Goudet, J. (2001). \emph{FSTAT, a program to estimate and test gene
diversities and fixation indices (version 2.9.3)} . Retrieved from
http://www.unil.ch/izea/softwares/fstat.html

Haffer, J. (2008). Hypotheses to explain the origin of species in
Amazonia. \emph{Brazilian Journal of Biology} , \emph{68} (4 SUPPL.),
917--947. doi: 10.1590/S1519-69842008000500003

Hakala, S. M., Seppä, P., \& Helanterä, H. (2019). Evolution of
dispersal in ants (Hymeno ptera: Formicidae): a review on the dispersal
strategies of sessile superorganisms Sanja. \emph{Myrmecological News}
,\emph{29} (March), 35--55. doi: 10.25849/myrmecol.news

Helmkampf, M., Gadau, J., \& Feldhaar, H. (2008). Population- and
sociogenetic structure of the leaf-cutter ant Atta colombica
(Formicidae, Myrmicinae). \emph{Insectes Sociaux} , \emph{55} (4),
434--442. doi: 10.1007/s00040-008-1024-3

Helms, J. (2018). The flight ecology of ants (Hymenoptera:
Formicidae).\emph{Myrmecological News} , \emph{26} , 19--30. doi:
10.25849/myrmecol.news

Hendry, A. P. (2004). Selection against migrants contributes to the
rapid evolution of ecologically dependent reproductive
isolation.\emph{Evolutionary Ecology Research} , \emph{6} (8),
1219--1236.

Hermelin, M. (2015). \emph{Landscapes and landforms of Colombia} (M.
Hermelin, Ed.). Springer International Publishing.

Hernández-Camacho, J. (1992). Caracterización geográfica de Colombia. In
G. Halffter (Ed.), \emph{La diversidad biológica de Iberoamérica
I}(Programa I, pp. 45--52). Acta Zoológica Mexicana.

Hölldobler, B., \& Wilson, E. O. (2011). \emph{The leafcutter ants:
civilization by instinct} (1st ed.). Nueva York, NY: W. W. Norton \&
Company.

Hoorn, C., Wesselingh, F. P., Ter Steege, H., Bermudez, M. A., Mora, A.,
Sevink, J., \ldots{} Antonelli, A. (2010). Amazonia through time: Andean
uplift, climate change, landscape evolution, and
biodiversity.\emph{Science} , \emph{330} , 927--931. doi:
10.1126/science.1194585

James, P. M. A., Coltman, D. W., Murray, B. W., Hamelin, R. C., \&
Sperling, F. A. H. (2011). Spatial genetic structure of a symbiotic
beetle-fungal system: Toward multi-taxa integrated landscape
genetics.\emph{PLoS ONE} , \emph{6} (10), 1--11. doi:
10.1371/journal.pone.0025359

Jombart, T., Devillard, S., \& Balloux, F. (2010). Discriminant analysis
of principal components: a new method for the analysis of genetically
structured populations. \emph{BMC Genetics} , \emph{11} (94), 1--15.
doi: 10.1371/journal.pcbi.1000455

Kattan, G. H., Franco, P., Rojas, V., \& Morales, G. (2004). Biological
diversification in a complex region: A spatial analysis of faunistic
diversity and biogeography of the Andes of Colombia. \emph{Journal of
Biogeography} , \emph{31} (11), 1829--1839. doi:
10.1111/j.1365-2699.2004.01109.x

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., \&
Mayrose, I. (2015). CLUMPAK: a program for identifying clustering modes
and packaging population structure inferences across K. \emph{Molecular
Ecology Resources} , \emph{15} (5), 1179--1191. doi:
10.1111/1755-0998.12387.CLUMPAK

Köppen, W. (1884). Die Wärmezonen der Erde, nach der Dauer der heissen,
gemässigten und kalten Zeit und nach der Wirkung der Wárme auf die
organische Welt betrachtet (The thermal zones of the earth according to
the duration of hot, moderate and cold periods and to the
impac.\emph{Meteorologische Zeitschrift} , \emph{1} , 215--226. doi:
10.1127/0941-2948/2011/105

Kronauer, D. J. C., Hölldobler, B., \& Gadau, J. (2004). Phylogenetics
of the new world honey ants (genus Myrmecocystus) estimated from
mitochondrial DNA sequences. \emph{Molecular Phylogenetics and
Evolution} , \emph{32} (1), 416--421. doi: 10.1016/j.ympev.2004.03.011

Lagomarsino, L. P., Condamine, F. L., Antonelli, A., Mulch, A., \&
Davis, C. C. (2016). The abiotic and biotic drivers of rapid
diversification in Andean bellflowers (Campanulaceae). \emph{New
Phytologist} , \emph{210} , 1430--1442. doi: 10.1111/nph.13920

Lee, C.-R., \& Mitchell-Olds, T. (2011). Quantifying effects of
environmental and geographical factors on patterns of genetic
differentiation. \emph{Molecular Ecology} , \emph{20} (22), 4631--4642.
doi: 10.5061/dryad.6rs51

Legendre, P., \& Fortin, M.-J. (2010). Comparison of the Mantel test and
alternative approaches for detecting complex multivariate relationships
in the spatial analysis of genetic data. \emph{Molecular Ecology
Resources} , \emph{10} (5), 831--844. doi:
10.1111/j.1755-0998.2010.02866.x

Librado, P., \& Rozas, J. (2009). DnaSP v5: A software for comprehensive
analysis of DNA polymorphism data. \emph{Bioinformatics} , \emph{25}
(11), 1451--1452. doi: 10.1093/bioinformatics/btp187

Luebert, F., \& Weigend, M. (2014). Phylogenetic insights into Andean
plant diversification. \emph{Frontiers in Ecology and Evolution}
,\emph{2} (27), 1--17. doi: 10.3389/fevo.2014.00027

Manel, S., \& Holderegger, R. (2013). Ten years of landscape
genetics.\emph{Trends in Ecology and Evolution} , \emph{28} (10),
614--621. doi: 10.1016/j.tree.2013.05.012

Manel, S., Schwartz, M. K., Luikart, G., \& Taberlet, P. (2003).
Landscape genetics: combining landscape ecology and population
genetics.\emph{Trends in Ecology and Evolution} , \emph{18} (4),
189--197. doi: 10.1016/S0169-5347(03)00008-9

Matute, D. R. (2013). The role of founder effects on the evolution of
reproductive isolation. \emph{Journal of Evolutionary Biology}
,\emph{26} (11), 2299--2311. doi: 10.1111/jeb.12246

Meirmans, P. G. (2006). Using the AMOVA framework to estimate a
standardized genetic differentiation measure. \emph{Evolution}
,\emph{60} (11), 2399--2402.

Meirmans, P. G. (2015). Seven common mistakes in population genetics and
how to avoid them. \emph{Molecular Ecology} , \emph{24} , 3223--3231.
Retrieved from http://dx.doi.org/10.1111/mec.13243

Meirmans, P. G., \& Hedrick, P. W. (2011). Assessing population
structure: FSTand related measures. \emph{Molecular Ecology Resources}
,\emph{11} (1), 5--18. doi: 10.1111/j.1755-0998.2010.02927.x

Montoya-Lerma, J., Giraldo-Echeverri, C., Armbrecht, I., Farji-Brener,
A., \& Calle, Z. (2012). Leaf-cutting ants revisited: Towards rational
management and control. \emph{International Journal of Pest Management}
,\emph{58} (3), 225--247. doi: 10.1080/09670874.2012.663946

Moser, J. C. (1967). Mating activities of Atta texana (Hymenoptera,
Fromicidae). \emph{Insectes Sociaux} , \emph{XIV} (3), 295--312.

Mueller, U. G., Ishak, H. D., Bruschi, S. M., Smith, C. C., Herman, J.
J., Solomon, S. E., \ldots{} Bacci, M. J. (2017). Biogeography of
mutualistic fungi cultivated by leafcutter ants. \emph{Molecular
Ecology} , (November), 1--17. doi: 10.1111/mec.14431

Muñoz-Valencia, V., Kähkönen, K., Montoya-Lerma, J., \& Díaz, F. (2020).
Characterization of a New Set of Microsatellite Markers Suggests
Polygyny and Polyandry in Atta cephalotes (Hymenoptera:
Formicidae).\emph{Journal of Economic Entomology} , \emph{113} (6),
3021--3027. doi: 10.1093/jee/toaa200

Muñoz-Valencia, V., Vélez-Matínez, G. A., Montoya-Lerma, J., \& Díaz, F.
(2021). Role of the Andean uplift as an asymmetrical barrier to gene
flow in the neotropical leaf-cutting ant Atta
cephalotes.\emph{Biotropica} , (October), 1--14. doi: 10.1111/btp.13050

Noguerales, V., Cordero, P. J., \& Ortego, J. (2016). Hierarchical
genetic structure shaped by topography in a narrow-endemic montane
grasshopper. \emph{BMC Evolutionary Biology} , \emph{16} (96), 1--15.
doi: 10.1186/s12862-016-0663-7

Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P.,
McGlinn, D., \ldots{} Wagner, H. (2019). Vegan: Community Ecology
Package. R Package Version 2.2-0. In \emph{Community ecology package} .
Retrieved from https://cran.r-project.org/web/packages/vegan/vegan.pdf

Parisod, C., Trippi, C., \& Galland, N. (2005). Genetic variability and
founder effect in the pitcher plant Sarracenia purpurea (Sarraceniaceae)
in populations introduced into Switzerland: from inbreeding to
invasion.\emph{Annals of Botany} , \emph{95} (2), 277--286. doi:
10.1093/aob/mci023

Peakall, R., \& Smouse, P. (2012). GenAlEx 6.5: Genetic analysis in
Excel. Population genetic software for teaching and
research.\emph{Bioinformatics} , \emph{1} , 6--8.

Peel, M. C., Finlayson, B. L., \& McMahon, T. A. (2007). Update world
map of the Köppen-Geiger climate classification. \emph{Hydrology and
Earth System Sciences} , \emph{11} , 1633--1644. doi: 10.1002/ppp.421

Pérez-Escobar, O. A., Gottschling, M., Chomicki, G., Condamine, F. L.,
Klitgård, B. B., Pansarin, E., \& Gerlach, G. (2017). Andean mountain
building did not preclude dispersal of lowland epiphytic orchids in the
Neotropics. \emph{Scientific Reports} , \emph{7} (1), 1--10. doi:
10.1038/s41598-017-04261-z

Pinheiro dos Reis, E., Fernandes Salomão, T. M., de Oliveira Campos, L.
A., \& Garcia Tavares, M. (2014). Genetic diversity and structure of
Atta robusta (Hymenoptera, Formicidae, Attini), an endangered species
endemic to the Restinga ecoregion. \emph{Genetics and Molecular Biology}
, \emph{37} (3), 581--586. doi: 10.1590/S1415-47572014000400015

Pritchard, J. K., Stephens, M., \& Donnelly, P. (2000). Inference of
population structure using multilocus genotype data. \emph{Genetics}
,\emph{155} (2), 945--959. doi: 10.1111/j.1471-8286.2007.01758.x

Ramos-Onsins, S. E., \& Rozas, J. (2002). Statistical properties of new
neutrality tests against population growth. \emph{Molecular Biology and
Evolution} , \emph{19} (12), 2092--2100. doi: Doi 10.1093/Molbev/Msl052

Rousset, F. (2008). GENEPOP'007: A complete re-implementation of the
GENEPOP software for Windows and Linux. \emph{Molecular Ecology
Resources} , \emph{8} (1), 103--106. doi:
10.1111/j.1471-8286.2007.01931.x

Rull, V. (2011). Neotropical biodiversity: timing and potential
drivers.\emph{Trends in Ecology and Evolution} , \emph{26} (10),
508--513. doi: 10.1016/j.tree.2011.05.011

Salgado-Roa, F. C., Pardo-Diaz, C., Lasso, E., Arias, C. F., Solferini,
V. N., \& Salazar, C. (2018). Gene flow and Andean uplift shape the
diversification of Gasteracantha cancriformis (Araneae: Araneidae) in
Northern South America. \emph{Ecology and Evolution} , \emph{8} (14),
7131--7142. doi: 10.1002/ece3.4237

Schowalter, T. D., \& Ring, D. R. (2017). Biology and management of the
Texas leafcutting ant (Hymenoptera: Formicidae). \emph{Journal of
Integrated Pest Management} , \emph{8} (1), 1--8. doi:
10.1093/jipm/pmx013

Sexton, J. P., Hangartner, S. B., \& Hoffmann, A. A. (2014). Genetic
isolation by environment or distance: which pattern of gene flow is most
common? \emph{Evolution} , \emph{68} (1), 1--15. doi: 10.1111/evo.12258

Shafer, A. B. A., \& Wolf, J. B. W. (2013). Widespread evidence for
incipient ecological speciation: a meta-analysis of
isolation-by-ecology. \emph{Ecology Letters} , \emph{16} (7), 940--950.
doi: 10.1111/ele.12120

Simon, C., Frati, F., Beckenbach, A., Crespi, B., Liu, H., \& Flook, P.
(1994). Evolution, weighting, and phylogenetic utility of mitochondrial
gene sequences and a compilation of conserved polymerase chain reaction
primers. \emph{Annals of the Entomological Society of America}
,\emph{87} (6), 651--701. doi: 10.1093/aesa/87.6.651

Siqueira, F. F. S., Ribeiro-Neto, J. D., Tabarelli, M., Andersen, A. N.,
Wirth, R., \& Leal, I. R. (2017). Leaf-cutting ant populations profit
from human disturbances in tropical dry forest in Brazil. \emph{Journal
of Tropical Ecology} , \emph{33} (5), 337--344. doi:
10.1017/S0266467417000311

Slatkin, M. (1993). Isolation by distance in equilibrium and
non-equilibrium populations. \emph{Evolution} , \emph{47} (1), 264--279.
Retrieved from http://www.jstor.org/stable/2410134

Smith, C. C., Weber, J. N., Mikheyev, A. S., Roces, F., Bollazzi, M.,
Kellner, K., \ldots{} Mueller, U. G. (2018). Landscape genomics of an
obligate mutualism: discordant population structures between a
leafcutter-ant and its fungal cultivars. \emph{BioRxiv} , 1--35. doi:
10.1101/458950

Sobel, J. M. (2014). Ecogeographic isolation and speciation in the genus
Mimulus. \emph{The American Naturalist} , \emph{184} (5), 565--579. doi:
10.1086/678235

Turchetto-Zolet, A., Pinheiro, F., Salgueiro, F., \& Palma-Silva, C.
(2013). Phylogeographical patterns shed light on evolutionary process in
South America. \emph{Molecular Ecology} , \emph{22} , 1193--1213. doi:
10.1111/mec.12164

VanOosterhout, C., Hutchinson, W. F., Wills, D. P. M., \& Shipley, P.
(2004). MICRO-CHECKER: Software for identifying and correcting
genotyping errors in microsatellite data. \emph{Molecular Ecology Notes}
, \emph{4} (3), 535--538. doi: 10.1111/j.1471-8286.2004.00684.x

Vasconcelos, H. L. (1997). Foraging activity of an Amazonian
leaf-cutting ant: Responses to changes in the availability of woody
plants and to previous plant damage. \emph{Oecologia} , \emph{112} (3),
370--378. doi: 10.1007/s004420050322

Wang, I. J., \& Bradburd, G. S. (2014). Isolation by
environment.\emph{Molecular Ecology} , \emph{23} (23), 5649--5662. doi:
10.1111/mec.12938

Wang, I. J., Glor, R. E., \& Losos, J. B. (2013). Quantifying the roles
of ecology and geography in spatial genetic divergence. \emph{Ecology
Letters} , \emph{16} (2), 175--182. doi: 10.1111/ele.12025

Wasko, A. P., Martins, C., Oliveira, C., \& Foresti, F. (2003).
Non-destructive genetic sampling in \selectlanguage{english}fish. An improved method for DNA
extraction from fish fins and scales. \emph{Hereditas} , \emph{138} ,
161--165.

Weber, J. N., Bradburd, G. S., Stuart, Y. E., Stutz, W. E., \& Bolnick,
D. I. (2016). Partitioning the effects of isolation by distance,
environment, and physical barriers on genomic divergence between
parapatric threespine stickleback. \emph{Evolution} ,\emph{71}
--\emph{2} , 342--356. doi: 10.1111/evo.13545

Wright, S. (1943). Isolation by distance. \emph{Genetics} , \emph{28} ,
114--138.

Zheng, C., Yang, F., Zeng, L., Vargo, E. L., \& Xu, Y. (2018). Genetic
diversity and colony structure of Tapinoma melanocephalum on the islands
and mainland of South China. \emph{Ecology and Evolution} , \emph{8}
(11), 5427--5440. doi: 10.1002/ece3.4065

\subsection*{DATA ACCESSIBILITY
STATEMENT}

{\label{data-accessibility-statement}}

The data that support the findings of this study will be openly
available in DataDryad after manuscript is accepted for publication.

\subsection*{AUTHORS' CONTRIBUTIONS}

{\label{authors-contributions}}

VMV participated in the conceptualization, ant sampling, data analysis,
original draft, reviewing, and editing of the manuscript. JML
participated in conceptualization, original draft, reviewing, and
editing of the manuscript. PS participated in data analysis, reviewing,
and editing of the manuscript. FD participated in conceptualization,
data analysis, original draft, reviewing, and editing of the manuscript.
All authors have read and approved the final version of the manuscript.

\subsection*{TABLES AND FIGURES}

{\label{tables-and-figures}}

\textbf{Table 1.} \textbf{Sampled} \textbf{populations of \emph{A.
cephalotes}.} Locations with their regional classifications (Pacific,
Andean 1, and Andean 2 regions from Colombia) and the number of nests
are indicated for each population.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
\textbf{Code} & \textbf{Location} & \textbf{Region} & \textbf{No. of
nests}\tabularnewline
\midrule
\endhead
BVT & Buenaventura & Pacific & 9\tabularnewline
QBD & Quibd\selectlanguage{ngerman}ó & Pacific & 15\tabularnewline
GUA & Guapi & Pacific & 20\tabularnewline
DAG & Dagua & Andean 1 & 10\tabularnewline
CAL & Cali & Andean 1 & 20\tabularnewline
PPY & Popay\selectlanguage{ngerman}án & Andean 2 & 14\tabularnewline
DOV & Dovio & Andean 2 & 19\tabularnewline
QBY & Quimbaya & Andean 2 & 16\tabularnewline
YOB & Yolomb\selectlanguage{ngerman}ó & Andean 2 & 19\tabularnewline
MNZ & Manizales & Andean 2 & 11\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2. Microsatellites and mitochondrial \emph{mtCOI} marker
variation for sampling locations of \emph{A. cephalotes} from the
Pacific, Andean 1, and Andean 2 regions of Colombia.} \emph{N} is the
sample size (number of nests). \emph{N\textsubscript{a}} : mean number
of alleles; \emph{AR} : mean allelic richness; \emph{H\textsubscript{e}}
: mean expected heterozygosity; \emph{F} \textsubscript{IS}: inbreeding
coefficient; \emph{HWE} : Hardy-Weinberg Equilibrium test; \emph{h} :
number of haplotypes; \emph{H\textsubscript{d}} : haplotype diversity;
\selectlanguage{greek}π: \selectlanguage{english}nucleotide diversity. (***) \emph{P} \textless{} 0.0001. Location
abbreviations are as presented in Table 1.\selectlanguage{english}
\begin{longtable}[]{@{}llllllllll@{}}
\toprule
& Nuclear data & Nuclear data & Nuclear data & Nuclear data & Nuclear
data & Mitochondrial data & Mitochondrial data & Mitochondrial data &
Mitochondrial data\tabularnewline
\midrule
\endhead
& N & N\textsubscript{a} & AR & H\textsubscript{e} &
\emph{F}\textsubscript{\textbf{IS}} & N & h & \textbf{\%
\emph{H\textsubscript{d}}} & \selectlanguage{greek}\textbf{\% π}\selectlanguage{english}\tabularnewline
\emph{Pacific} & & & & & & & & &\tabularnewline
All pop n = 3 & 44 & 6.95 & 6.56 & 0.57 & 0.01 & 42 & 5 & 67 &
0.4\tabularnewline
BVT & 9 & 6.08 & 6.02 & 0.59 & -0.06 & 9 & 2 & 22 & 0.2\tabularnewline
QBD & 15 & 7.46 & 6.97 & 0.58 & 0.00 & 13 & 4 & 60 & 0.4\tabularnewline
GUA & 20 & 7.31 & 6.68 & 0.55 & 0.04 & 20 & 4 & 56 & 0.4\tabularnewline
\emph{Andean 1} & & & & & & & & &\tabularnewline
All pop n = 2 & 30 & 4.62 & 4.38 & 0.49 & -0.01 & 29 & 6 & 82 &
0.5\tabularnewline
DAG & 10 & 3.92 & 3.85 & 0.46 & -0.01 & 10 & 3 & 51 & 0.4\tabularnewline
CAL & 20 & 5.31 & 4.91 & 0.51 & -0.01 & 19 & 4 & 76 & 0.3\tabularnewline
\emph{Andean 2} & & & & & & & & &\tabularnewline
All pop n = 5 & 79 & 4.74 & 4.51 & 0.53 & -0.01 & 75 & 6 & 53 &
0.2\tabularnewline
PPY & 14 & 3.69 & 3.58 & 0.54 & -0.01 & 14 & 3 & 62 & 0.2\tabularnewline
DOV & 19 & 4.38 & 4.18 & 0.48 & -0.05 & 19 & 3 & 29 & 0.1\tabularnewline
QBY & 16 & 4.77 & 4.58 & 0.52 & 0.02 & 15 & 2 & 25 & 0.1\tabularnewline
YOB & 19 & 6.38 & 5.83 & 0.61 & -0.01 & 16 & 1 & 0 & 0.0\tabularnewline
MNZ & 11 & 4.46 & 4.37 & 0.50 & 0.02 & 11 & 3 & 35 & 0.3\tabularnewline
\textbf{Total n = 10} & 153 & 10.38 & 7.06 & 0.60 & -0.01 & 146 & 10 &
66 & 0.3\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 3.} \textbf{Hierarchical AMOVA for microsatellites
and\emph{mtCOI} markers.} The analyses included two intermediate levels
of variation following: 1) isolation by barrier (IBB), with populations
separated by the western range (Pacific and Andean), and 2) isolation by
environment (IBE), with populations classified into three regions
following environmental conditions (Pacific, Andean 1, and Andean 2).
The proportion of variation (\%) and
associated\emph{F\textsubscript{ST}} and \selectlanguage{greek}\emph{Φ} \selectlanguage{english}are given for each
hierarchical level.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllllll@{}}
\toprule
\textbf{Model} & & & & Microsatellites & Microsatellites &
Microsatellites & Microsatellites & mtCOI & mtCOI & mtCOI\tabularnewline
\midrule
\endhead
& & & \textbf{~} \emph{F}\textsubscript{ST} & \textbf{~}
\emph{F}\textsubscript{ST} & \textbf{~} \emph{F}\textsubscript{ST} &
\textbf{~} \emph{F}\textsubscript{ST} & \selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{ST} &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{ST} & \selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{ST} &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{ST}\tabularnewline
& \textbf{Source of variation} & \textbf{Source of variation} &
\textbf{df} & \textbf{df} & \textbf{\%} & F\textsubscript{ST} &
\textbf{df} & \textbf{df} & \textbf{\%} & \selectlanguage{greek}\emph{Φ}\selectlanguage{english}\tabularnewline
\textbf{1)} IBB: Pacific populations and Andean populations & Regions &
Regions & 1 & 1 & 2 & \emph{F}\textsubscript{CT}: 0.02* & 1 & 1 & 0 &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{CT}: -0.08\tabularnewline
& Pops (regions) & Pops (regions) & 8 & 8 & 9 &
\emph{F}\textsubscript{SC}: 0.09* & 8 & 8 & 33 &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{SC}:0.33*\tabularnewline
& Within pops & Within pops & 296 & 296 & 90 &
\emph{F}\textsubscript{ST}:0.10* & 136 & 136 & 67 &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{ST}: 0.28*\tabularnewline
& Total & Total & 305 & 305 & 100 & & 145 & 145 & 100 &\tabularnewline
\textbf{2)} IBE: Pacific, Andean 1, and Andean 2 populations & Regions &
Regions & 2 & 2 & 4 & \emph{F}\textsubscript{CT}: 0.04* & 2 & 2 & 0 &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{CT}: -0.06\tabularnewline
& Pops (regions) & Pops (regions) & 7 & 7 & 7 &
\emph{F}\textsubscript{SC}: 0.07* & 7 & 7 & 33 &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{SC}:0.34*\tabularnewline
& Within pops & Within pops & 296 & 296 & 89 &
\emph{F}\textsubscript{ST}: 0.11* & 136 & 136 & 67 &
\selectlanguage{greek}\emph{Φ}\selectlanguage{english}\textsubscript{ST}: 0.30*\tabularnewline
& Total & Total & 305 & 305 & 100 & & 145 & 145 & 100 &\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 4. Results of model selection as tested through the RDA
analysis under IBD, IBE, and IBB scenarios of genetic divergence.}\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
\textbf{MODEL} & \textbf{F} & P & \textbf{AIC}\tabularnewline
\midrule
\endhead
One model & & &\tabularnewline
IBD model & 1.17 & 0.30 & 1.33\tabularnewline
IBD model selected (PCNM1) & 2.71 & 0.01* & 0.49\tabularnewline
IBE model (T + H + P + E + C) & 1.46 & 0.07 & 1.65\tabularnewline
IBE model selected (T + P) & 1.81 & 0.03* & 1.23\tabularnewline
IBB model (Andes) & 1.31 & 0.22 & 1.89\tabularnewline
Two models combined & & &\tabularnewline
IBE-IBB (T + P + Andes) & 1.69 & 0.03* & 1.28\tabularnewline
IBE-IBD (T + P + PCNM1) & 2.16 & [?] 0.01* & 0.09\tabularnewline
IBB-IBD (Andes + PCNM1) & 2.17 & [?] 0.01* & 0.58\tabularnewline
Three models combined & & &\tabularnewline
Whole model (IBD + IBE + IBB) & 2.17 & [?] 0.01* & -0.66\tabularnewline
Conditioned models by IBD & & &\tabularnewline
IBE model (T + P) \textasciitilde{} PCNM1 & 1.65 & 0.04* &
0.09\tabularnewline
IBB model (Andes) \textasciitilde{} PCNM1 & 1.48 & 0.13 &
0.58\tabularnewline
Whole model (IBD + IBE + IBB) \textasciitilde{} PCNM1 & 1.74 & 0.02* &
-0.66\tabularnewline
\bottomrule
\end{longtable}

T: Temperature; H: Humidity; P: Precipitation; E: Elevation; C: Climate.
(*) Indicates significant value, \emph{P} [?] 0.01. AIC null model: 1.41.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1. Study area and sampling locations of \emph{Atta
cephalotes}} \textbf{in Colombia.} Different symbols and colors indicate
sampling locations following regional classification: Blue triangles for
sampling locations in the Pacific region (BVT, QBD, and GUA), violet
squares for the Andean 1 region (DAG and CAL), and orange circles for
the Andean 2 region (PPY, DOV, QBY, YOB, MNZ) (Table 1, and
Supplementary table ST1).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2.} \textbf{Patterns of genetic diversity and haplotype
network for} \emph{A. cephalotes} \textbf{sampled in three different
regions of Colombia}. \textbf{A.} Genetic diversity parameters for
microsatellite markers (N\textsubscript{a}: mean number of different
alleles; A\textsubscript{R}: mean allelic richness; H\textsubscript{E}:
mean expected heterozygosity). \textbf{B.} Genetic diversity parameters
for the \emph{mtCOI} gene (h: number of haplotypes; \selectlanguage{greek}π: \selectlanguage{english}nucleotide
diversity; H\textsubscript{d}: haplotype diversity).
\textbf{C.}Haplotype network based on 146 sequences of the \emph{A.
cephalotesmtCOI} gene. Numbers on the connecting lines denote the number
of mutations.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3. Heatmap for pairwise
\emph{F}\textsubscript{ST}comparisons of \emph{Atta cephalotes}
populations. A.}\emph{F} \textsubscript{ST} comparisons based on nuclear
microsatellites and, \textbf{B.} \selectlanguage{greek}\emph{Φ} \selectlanguage{english}\textsubscript{ST} comparisons
based on\emph{mtCOI} . Red squares denote intra-region comparisons.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 4.} \textbf{Population structure of \emph{Atta
cephalotes} as obtained from the STRUCTURE and DAPC analyses.
A.}Discriminant analysis of principal components (DAPC) scatterplot
between samples. The key describes the colors attributed to each region
and inertia ellipses describe the general distribution of points.
Eigenvalues for each PC axis are shown (PC1, vertical; PC2, horizontal).
The number of PCA axes retained in each DAPC analysis is shown in the
bottom-right inset (gray bars). \textbf{B.} Multilocus genotype
clustering of \emph{A. cephalotes} populations, using STRUCTURE and
following Evanno et al. (2005). Each region is characterized by color
and each location is represented by a vertical bar. Locations are
organized by region, starting with Pacific (BVT -- GUA) followed by
Andean 1 (DAG and CAL), and Andean 2 (PPY -- MNZ). A different
simulation was run for the southern cluster. \textbf{C.} Map for
sampling localities according to the STRUCTURE clustering: Pacific (BVT,
QBD, and GUA in blue), Andean south (DAG, CAL, and PPY in violet), and
Andean north (DOV, QBY, YOB, and MNZ in orange).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5.} \textbf{Variation partitioning and model comparison
of isolating mechanisms driving nuclear genetic differentiation
(\emph{F}\textsubscript{ST}) in \emph{Atta cephalotes}.} Variation
partitioning analysis based on RDA results of the full model (IBE + IBB
+ IBD) into environmental (IBE: temperature and precipitation), barrier
(IBB: Andes mountains), and spatial (IBD: geographic distance)
components. Each circle represents the variation explained for each
mechanism, and their overlap is the fraction of shared variation. Model
selection and significance are indicated by the AIC and \emph{P}-values.
The IBE model remained significant after controlling for IBD
(Conditioned model: F: 1.66, \emph{P}: 0.04, AIC: 0.09), while IBB was
not significant when accounting for IBD (Conditioned model: F:
1.49,\emph{P}: 0.11, AIC: 0.56).

\selectlanguage{english}
\FloatBarrier
\end{document}

