Animal collection and culturing

We discovered the M. emersoni wingless phenotype in the Chiricahuas mountain range in 1999. We found that this wingless phenotype exists on four other Sky Island mountains in Southeastern Arizona, USA and we collected colonies from 14 sites (Additional file 1: Table S3). We preserved colonies in 95 % ethanol or kept them alive in plastic boxes with glass test tubes filled with water constrained by a cotton plug. We fed colonies with a mix of crickets, mealworms and Bhatkar-Whitcomb diet [105]. The colonies were maintained at 27 ° C and 70 % humidity with a day/night cycle of 12 hours. We induced production of winged and wingless queens by exposing them to a cycle of temperature range: the temperature of the growth chamber was reduced by 1 ° C for 17 days until is reached 10 ° C. We kept the chamber at 10 ° C for 2 weeks, and then increased the temperature gradually by 1 ° C per day for 17 days until the original conditions were restored.

Statistical analyses

We performed all statistical analyses using R, unless otherwise stated.

Inferring modern and historical population demography

Sequencing and genotyping

We extracted DNA from 318 M.emersoni queens, each from a different colony, using thoracic tissue from preserved individuals in ethanol. We amplified by PCR two mitochondrial genes (cytochrome oxidase 1 (CO1) and 2(CO2)), and sequenced them in both directions on an ABI3730 (for primers, see [106]) for a total of 1350 base pairs (bp) for 70 individuals. We also amplified and sequenced 1020 bp from a nuclear gene, hymenoptaecin, including one intron. For each Sky Island, we isolated alleles from hymenoptaecin to create a pool of DNA for each Sky Island by mixing DNA from about 20 queens from different colonies. This method has been shown to be accurate to estimate allele frequencies of a population [107, 108]. We amplified all alleles in the DNA pool simultaneously by PCR and cloned them into a PGEM-T vector. We picked 15 bacterial clones per mountain and we sequenced them with M13 primers in both directions. This allowed us to obtain the exact haplotype phase for every nuclear allele sampled and sequenced. All sequences were assembled and aligned with Geneious Pro V.5.4.3 and checked by eye. GenBank accession numbers are KT356282-KT356545.

We also genotyped all individuals for Amplified Fragment Length Polymorphisms (AFLP) markers, as described in [109]. We used 4 primer combinations to generate 157 unambiguous loci: M CTA and E ACC , M CTA and E ACG , M CAG and E ACC , and M CAG and E ACG . We ran the PCR products on a LiCor DNA analyser 4300 and the same researcher scored them blindly using the accompanying Saga software. To estimate the error rate in the procedure, about 10 % of the samples were amplified and scored twice, and we recovered a similar error rate to what has been reported in other studies [110].

Estimation of modern gene flow

We estimated F ST values with the Bayesian method developed by [111] and implemented in AFLP-surv [112]. We removed from the analysis sites where we collected less than 15 individual queens from different colonies. P-values for the significance of each F ST were estimated after applying Boneferroni corrections for multiple comparisons. Pairwise F ST estimates for all 14 sites are reported in Additional file 1: Table S5. Stars represent parameters that could not be estimated with confidence. To test for the presence of recent gene flow, we performed a Mantel test [113] to calculate the correlation between the three-dimensional geographic distance and the genetic differentiation (F ST estimates) for all pairwise comparisons among sites. We expected a significant correlation to occur if restricted but non-zero gene flow existed among Sky Islands, or if Sky Islands diverged recently from a continuous ancestral population without recent gene flow. On the other hand, we expected a non-significant correlation for the two following demographic scenarios: a complete panmixia among most sites, or a highly reduced gene flow for a long time among most sites. We tested the Mantel coefficient for significance using 10,000 random permutations on the matrices.

Estimation of historical gene flow and times of divergence

We used the Bayesian search strategy implemented in MIGRATE-n 3.2.15 [114] to estimate the effective population sizes, the average historical migration rates, the historical migration rates and the coalescence times for each population. We also evaluated 2 models of population subdivision by comparing their marginal likelihoods using the modified thermodynamic integration (TH) implemented in MIGRATE-n [115]. The first model considered all sampling sites as part of one panmictic unit while the second model considered each mountain range as independent units (limited or absence of gene flow). Each MIGRATE-n run consisted of 10 replicates of 4 MCMC heated short chains (heating scheme: 1.0, 1.2, 3.0, 6.0) and one long chain of 5,000,000 sampled trees where 500,000 burn-in trees were discarded to ensure proper sampling of the parameter space. All parameter estimates showed narrow posterior probability distributions, indicating high confidence in the estimates. All comparisons involving the Chiricahuas as the sending population show clear outlier values and were thus not considered. Θ values from the haploid and maternally-transmitted mitochondrial loci were corrected for comparison with a diploid locus to obtain an estimate of 4Ne μ. We used the direct estimate of mtDNA mutation rate calculated by [116] and the immune gene mutation rate estimate in Drosophila from [117] to calculate migration rates and times of divergence. MIGRATE-n also implements the Skyline plot method for detecting changes in migration rates over time, thus allowing the reconstruction of historical migrant exchanges between mountain ranges. The analysis was run separately for each pair of geographically adjacent mountain ranges (Chiricahuas-Huachucas, Huachucas-Catalinas, Catalinas-Pinals, Pinals-Pinaleños and Pinaleños-Chiricahuas) with the same run parameters as described above.

Testing for competing evolutionary scenarios of population splits

We used Approximate Bayesian Computation implemented in the program DIYABC v.1.0.4.45 beta [118] to test competing hypotheses of lineage divergence. Instead of estimating the likelihood of a model from an MCMC approach, DIYABC uses approximate Bayesian computation [119], where similarity of summary statistics between observed and simulated data sets are compared. Data can be simulated for several scenarios, where population sizes, times of divergence, population size change and admixture can be specifically defined.

In order to focus on testing the most probable scenarios, we built up our hypotheses based on the results of the Skyline plot analysis from MIGRATE-n (Additional file 1: Table S4 and Additional file 1: Figure S1). We compared three divergence scenarios; (A) a simultaneous isolation of the five mountain ranges from a common ancestor, (B) a sequential divergence of Sky Islands based on the results of the Skyline plot analysis using as priors the estimated times where no more gene flow is observed between two Sky Islands, and (C) a third scenario based on scenario B, but with the Huachucas and the Catalinas interchanged on the topology, thereby affecting the relationships within the Southern and Northern Sky Islands (Additional file 1: Figure S2).

The prior probability uniform distribution on the oldest divergence event (t4) ranged from 10,000 to 200,000 generations in the past and the prior probability uniform distribution on the effective population (Ne) size ranged from 1000 to 200,000 (coherent with MIGRATE-n estimates of Ne). The sequence of subsequent splits was forced to be kept in this specific order when simulating the data and the time interval prior probability distribution for each of these respective splits were t4 =1 to 50,000, t3 =10,000 to 70,000, t2 =25,000 to 100,000 and t1 =10,000 to 200,000 generations in the past. The prior on Ne was kept constant through time and was sampled from a uniform distribution ranging from 1000–100,000. Most importantly, three conditions were imposed for the sampling of parameters for the population splitting to occur in the specified order: t1 ≥ t2, t2 ≥ t3 and t3 ≥ t4. 300,000 simulations were performed and the direct estimate method was used to estimate the posterior probabilities of each scenario.

Parallel evolution of M.emersoni populations

Inferring mitochondrial haplotype network

To get insight about whether the wingless phenotype evolved repeatedly in parallel or has a single origin, we first infered genealogical realtionships among queen haplotypes across the Sky Islands populations. A mitochondrial DNA network was constructed using the R package ape. Divergent haplotypes (12 and 20 steps away from the central haplotype were discovered in the Chiricahuas but left out of the analysis for simplicity of the representation of the network.

Neighbor-joining tree

A neighbor-joining tree was constructed using AFLP data with the phylogenetic package Phylip [120]. The tree was constructed using sequentially the executables Genedist, Neighbor, Seqboot, Consense and Drawtree. Bootstrap values over 70 are reported on the tree.

Redundancy analyses

We first calculated the pairwise mtDNA genetic distance between each pairs of queens and performed a multidimensional scaling to obtain a set of synthetic variables that best represented the pairwise distances between records. To quantify habitat type in a discrete manner, we performed a principal component analysis using all environmental variables on all collection sites, assigning a habitat type (1 to 5) to sites with high similarity. We built three explanatory matrices: collection site, phenotype (winged or wingless) and environment (habitat type). We partitioned the variation in genetic distance with respect to the three explanatory matrices using a redundancy analysis ordination.

Adaptation to climatic changes

Tests for temperature tolerance

Heat resistance of worker ants was measured using a knockdown assay. About 20 individuals workers were collected from inside the nest of each colony, placed into 8 mL glass vials, and incubated at 50 ° C for one hour into a hybridization oven. Resistance was scored for each ant as the time taken to be knock down and mortality after treatment was recorded for each sample. Cold resistance was evaluated as recovery from a chill coma. About 20 individuals workers were collected from inside the nest of each colony, placed into small petri dishes, and immersed for 16 hours into a box filled with ice and water. This treatment led to rapid immobilization of all individuals. The temperature was kept constant around 4 ° C. Petri dishes were then exposed to ambient temperature and recovery was scored for each worker as the time taken for muscular coordination to come back.

Environmental factors the winged and wingless phenotypes may be adapting to along the ecological gradients

Explanatory variables for queen phenotype distribution

We calculated the percentage of wingless queens present at each site. This value represents the proportion of wingless queens that were successful in founding or persisting in a colony, and thus, have survived the dispersal/founding selection period. For each site, we measured the latitude and elevation in the field, and we obtained climatic variables from WorldClim Atlas and the vegetational growth average and vegetational growth peak were obtained from the National Atlas (Additional file 1: Table S1). The Worldclim Atlas and the National Atlas data have a 1 km2 resolution, which is, from our observations, of the same magnitude as the area of our average sampling sites.

We performed an information-theoretic approach to evaluate the performance of 12 models explaining the frequency distribution of winged and wingless queens along the ecological gradients. This approach, widely used in ecology, does not reject or accept models, but rather ranks a priori models and provides posterior probabilities that allow for the interpretation of the performance of each model. Using the findings of previous studies in the literature, we constructed 10 biologically plausible models that could explain which ecological factor(s) have the most influence on distribution of wingless queens along the ecological gradients. We used the method and R code of [121] to fit and rank models: we fitted models using the R function lm, and we calculated the sample-size corrected Akaike Information Criterion (AICc) for each model. We controlled for population differentiation using the score of the first component from a principal component analysis on the mtDNA pairwise genetic distances among all individual queens. We ranked the models based on AICc and calculated their posterior and conditional probabilities for further interpretation.

Interruptions of the gene network controlling wing development

Sample collection

We collected queen larvae at late last larval instar, just prior to the prepupal stage, following the criteria of [122]. We fixed larvae for 2h in a 4 % formaldehyde solution in PEM buffer [123], and dissected them under a Zeiss Discovery V12 Stereomicroscope to expose their thoracic imaginal discs by removing obstructive tissues.

Wing imaginal disc morphology and cell division

We calculated surface area of the forewing imaginal disc and of a leg imaginal disc (in μm2) using AxioVision software (Carl Zeiss Canada Ltd., Toronto, Ontario, Canada). For all samples that were undamaged we counted the number of Phosphohistone expressing cells in the forewing imaginal disc from a z-stack image. To test for a difference in static allometry of leg and wing disc among populations, we performed a Standardized Major Axis bivariate line-fitting analysis (SAM) and tested for slope differences [124] (Additional file 1: Figure S6A). We also performed an analysis of variance on the number of cells undergoing mitosis per wing disc unit area for each population (Fig. 4).

Expression patterns of genes of the wing patterningnetwork

Immunohistochemistry: we performed antibody stainings for Engrailed (En), Ultrabithorax (Ubx), Cut (Cut) and Extradenticle (Exd) proteins following the protocol of [123]. Antibody stainings for Myosin enhancing factor-2 (Mef2) and Phosphohistone-3 (PH3) were conducted similarly, except that the tissues were treated with PAT 1 % instead of PTW 1 % to increase the antibody penetration in the tissue. We systematically included winged queen larvae in each staining reaction containing wingless queen larvae as positive controls and we over developed the staining reaction for larvae of both winged and wingless queens to confirm that the absence of expression was real and not an artifact. We imaged the stained tissues with a Zeiss microscope using the AxioVision software. Expression of genes in the network that controls wing development is highly conserved across winged castes of ants [63, 65]. Nonetheless, we confirm that winged queens from several collection sites do not show any expression differences (data not shown).

Whole mount in-situ hybridization: We isolated a fragment of the gene wingless by PCR using the primers of [63] on cDNA synthesized from embryonic and larval RNA. We cloned and sequenced this fragment to confirm its identity. Genbank accession number is KT361189. We followed the whole mount in-situ protocol from [125].

3D X-ray microtomography

In order to describe in more detail the differences in external and internal anatomy of the thorax, we made 3D X-ray microtomography images (microCT scans) of winged queens and 11 wingless queens of different geographic origins. Ants were stained with elemental iodine (1 % w/v in absolute ethanol) to enhance the X-ray contrast of the soft tissues [126]. The scans were performed on an Xradia MicroCT 3D imaging system with voxel sizes between 1.8 and 4.2 μm. We examined sagittal, coronal, and transverse sections for the presence of direct and indirect flight muscles. We applied a surface rendering algorithm using OsiriX to visualize the 3D reconstruction of the external anatomy of the thorax.