Skip to main content
  • Research Article
  • Open access
  • Published:

Past climate change on Sky Islands drives novelty in a core developmental gene network and its phenotype



A fundamental and enduring problem in evolutionary biology is to understand how populations differentiate in the wild, yet little is known about what role organismal development plays in this process. Organismal development integrates environmental inputs with the action of gene regulatory networks to generate the phenotype. Core developmental gene networks have been highly conserved for millions of years across all animals, and therefore, organismal development may bias variation available for selection to work on. Biased variation may facilitate repeatable phenotypic responses when exposed to similar environmental inputs and ecological changes. To gain a more complete understanding of population differentiation in the wild, we integrated evolutionary developmental biology with population genetics, morphology, paleoecology and ecology. This integration was made possible by studying how populations of the ant species Monomorium emersoni respond to climatic and ecological changes across five ‘Sky Islands’ in Arizona, which are mountain ranges separated by vast ‘seas’ of desert. Sky Islands represent a replicated natural experiment allowing us to determine how repeatable is the response of M. emersoni populations to climate and ecological changes at the phenotypic, developmental, and gene network levels.


We show that a core developmental gene network and its phenotype has kept pace with ecological and climate change on each Sky Island over the last 90,000 years before present (BP). This response has produced two types of evolutionary change within an ant species: one type is unpredictable and contingent on the pattern of isolation of Sky lsland populations by climate warming, resulting in slight changes in gene expression, organ growth, and morphology. The other type is predictable and deterministic, resulting in the repeated evolution of a novel wingless queen phenotype and its underlying gene network in response to habitat changes induced by climate warming.


Our findings reveal dynamics of developmental gene network evolution in wild populations. This holds important implications: (1) for understanding how phenotypic novelty is generated in the wild; (2) for providing a possible bridge between micro- and macroevolution; and (3) for understanding how development mediates the response of organisms to past, and potentially, future climate change.


How populations differentiate within species is a fundamental problem in evolutionary biology that is key for uncovering the processes that generate biological diversity, including microevolution, speciation, and the emergence of phenotypic novelty [16]. Traditionally, biologists have used the tools of population and quantitative genetics, two fields central to the modern evolutionary synthesis, to study population differentiation [1, 79]. These fields have significantly advanced our knowledge of variation and change in the frequency of alleles and phenotypes (quantitative and discrete) within and between populations [10]. They have also uncovered signatures of natural selection and genetic drift and are being used to identify loci responsible for adaptive and non-adaptive phenotypes driving the evolution of populations [8, 1114]. However, a largely unexplored dimension of this problem is what role, if any, does organismal development play in the process of population differentiation in the wild [3, 15, 16].

Organismal development integrates the action of gene regulatory networks with environmental inputs to generate the phenotype [3, 1719], and therefore, may play a key role in facilitating phenotypic differentiation of populations and species exposed to ecological changes [3, 17, 19]. Major advances in the field of evolutionary developmental biology, such as the discovery of Hox genes, have revealed an unexpected degree of evolutionary conservation of developmental regulatory genes across the animal kingdom [7, 20]. These regulatory genes, which are transcription factors and signalling molecules, are assembled in hierarchically organized networks [18, 20]. A primary function of developmental networks is body plan formation [18, 20]. For example, the expression, structure, function, and regulation of the developmental genes Sonic hedgehog and patched have been conserved for hundreds of millions of years across vertebrates and invertebrates [21]. Both Sonic hedgehog and its downstream target patched regulate anterior/posterior patterning in the developing limb of both fruit flies and chickens [21]. The high degree of conservation and structure of developmental networks suggests that these networks may bias the variation that selection can act upon. Such bias may facilitate repeatable phenotypic responses when populations are exposed to similar environmental inputs or encounter similar ecological changes [1, 3, 22, 23]. Here we use an integrative approach that combines multiple levels of organization (gene networks, development, and phenotype) as well as multiple fields (population genetics, paleoecology and ecology, morphology, and evolutionary developmental biology) to gain a more complete understanding of population differentiation in the wild.

Our ability to integrate all of these levels and approaches in a single study was made possible by using recently diverged populations within the ant species Monomorium emersoni along five ‘Sky Islands’ in Arizona. Sky Islands are a group of mountain ranges in the American Southwest across the Arizona-Mexico border that are isolated by large areas of deserts with limited genetic exchange between them [24, 25]. From the perspective of the evolutionary biologist, Sky Islands represent a replicated natural experiment, in which each Sky Island is a natural laboratory of high evolutionary potential. Isolation of the Sky Islands has facilitated the evolution of several endemic species: several squirrels, plants, ants, more than 60 species of land snails, as well as many endangered species that are found only in the Sky Islands of Arizona [2629]. In essence, this archipelego of Sky Islands is comparable to oceanic islands isolated by sea water, like those in the Galapagos, Hawaii and the Caribbean islands [27] (Fig. 1 a and b).

Fig. 1
figure 1

The Arizona Sky Islands. a Satellite image of North America showing the location of the Arizona Sky Islands (yellow box). b A picture taken on a Sky Island from high elevation. Across the desert at the horizon, a second Sky Island can be distinguished. Alternative phenotypes of Monomorium emersoni queens: c winged queen and d wingless queen. e Map of Southeastern Arizona showing the relative location of the five Sky Islands indicating our sampling sites (black dots) and 500m topographic contour lines. f Representation of the highly conserved gene network that controls wing development in Drosophila melanogaster. Arrows indicate activation and bars indicate repression. Genes examined in this study are shown in black

The paleoecological record of the Sky Islands is one of the most complete records available [24, 25]. Numerous studies have reconstructed the major ecological and climatic changes that have occurred in the region by using the stratification of pollen grains and middens of pack rats in the fossil record. These studies show that during the Pleistocene glacial periods 90,000–20,000 years BP when atmospheric temperature around the region was substantially cooler, there was a largely continuous forest landscape connecting Sky Island mountain ranges [24, 25]. The region then began warming up between 20,000–10,000 years BP, and atmospheric temperatures increased substantially, resulting in the formation of vast deserts that isolated the Sky Islands [24, 25]. Species that occurred at lower elevations in the intervening valleys during the Pleistocene glacial periods, such as Ponderosa pine trees and Douglas fir, are currently only found at high-elevations. These past climatic changes in the region have generally caused shifts in elevation, expansions, and contractions of populations in numerous species [2832]. Each Sky Island is now considered to be one of the steepest terrestrial ecological gradients in North America, from low-elevation deserts, through mid-elevation oak-juniper woodlands, to high elevation coniferous forests [27]. Several ecological factors vary steeply along these altitudinal clines, in particular temperature, precipitation and habitat fragmentation [3335]. Environmental changes that would normally be observed across a latitudinal gradient stretching thousands of kilometres can be observed across a few hundred meters of elevational change in the Sky Islands [26, 27, 36]. This poses a unique adaptive challenge for organisms on Sky Islands that have distributional ranges that encompass more than one altitudinal zone [26, 27, 36]. These steep gradients therefore behave as powerful sensors of climatic and environmental changes because noticeable changes in communities occur rapidly and over short geographical distances [3335].

Ants are an ecologically dominant group [37] that have long been used as bioindicators of environmental change [38, 39] and climate warming [40, 41]. In most ant species, the colony is made up of a single winged queen that performs most of the reproduction in the colony and her wingless workers that perform almost all the other tasks [37]. Throughout most of the year, the queen produces eggs that develop into wingless workers, but in response to specific environmental conditions associated with spring time, such as the gradual increase in daylight and temperature, they begin producing virigin queens and males that use their fully functional wings to disperse and participate in mating flights away from the mother colony. Immediately after these mating flights, males die while queens tear off their wings to found a new nest underground [37]. This life history strategy is known as ‘independent colony foundation’ [4245]. In some ant species, however, an alternative life history strategy called ‘dependent colony foundation’ has evolved [4245]. In these species, there are multiple queens in a single ant colony in which some or all of the queens are wingless. Virgin queens that are wingless mate on the ground and disperse on foot aided and accompanied by wingless workers from the mother colony. This strategy dramatically decreases mortality of wingless queens, but at the expense of long distance dispersal. Wingless queens have evolved multiple times independently within ants [4345], suggesting that wingless queens represent an adaptive life-history strategy. It takes less energy to produce wingless queens, and selection may favor the evolution of wingless queens in particular ecological conditions where aerial dispersal is risky [29, 46, 47].

Colonies of the ant species M. emersoni are distributed throughout the Sky Islands, where we discovered the presence of two alternative queen phenotypes, winged and wingless (Fig. 1 c and d) [48]. A single M. emersoni colony can be headed by one or multiple queens, ranging from 1 to more than 70, and colonies can be composed either exclusively by winged queens, exclusively by wingless queens, or a mixture of winged and wingless queens [48]. Current evidence suggests that determination of alternative queen phenotypes in a colony occurs during development and can be influenced by both environmental and genetic factors [43, 49, 50]. On each Sky Island, M. emersoni occurs sparsely at lower elevations, where its abundance is very low and is associated with creek or woody patches near mountain ranges [48]. At mid-elevations of each Sky Island, from 1500–2300 meters, M. emersoni occurs at high densities in oak-juniper woodlands, whereas at the highest elevations of each Sky Island, from 2300 m and above, M. emersoni occurs in patches of high density in coniferous forests dominated by Pinyon pine and Douglas fir [48].

Studying how population differentiation occurs in M. emersoni along several naturally replicated ecological gradients on five Sky Islands in Arizona provides a unique opportunity to ask whether organismal development facilitates repeatable, and thus, predictable responses to similar ecological and climate changes on each Sky Island. Furthermore, the interaction between ecology and organismal development in the context of Sky Islands is key for understanding population dynamics of gene networks and its impacts on phenotypic evolution on microevolutionary timescales, as well as for understanding how past changes in developmental gene networks may affect the future course of evolution in response to climate change.


We first used heat and cold shock experiments to determine whether colonies of M. emersoni have responded to the establishment of warmer conditions at the lower limit of their elevational range across all 5 Sky Islands during the last 10,000 years BP. We found that colonies from warmer sites show increased tolerance to high temperature, whereas colonies from colder sites show increased tolerance to cold temperature (Fig. 2 a and b). This suggests that individuals in colonies of M. emersoni have adapted to local temperature conditions along the ecological gradients. To determine how populations have changed their distribution over the last 100,000 years BP (Fig. 1 e, Additional file 1: Table S3), we reconstructed historical levels of gene flow using mitochondrial (mtDNA, CO1 and CO2) and nuclear (hymenotpaecin) sequences (Additional file 1: Table S4 and Additional file 1: Figure S1). Our coalescent (Additional file 1: Figure S1) and evolutionary (Additional file 1: Figure S2) simulations show a sequential model of isolation between adjacent Sky Islands. A first split between southern (Huachucas and Chiricahuas) and northern (Pinaleños, Catalinas, and Pinals, Fig 1 e and 3 a) Sky Islands occurred around 80,000 years BP and was followed by the complete isolation of the 5 mountain ranges between 20,000 and 10,000 years BP (Additional file 1: Table S5, Figure S1 and S2). These results indicate that some Sky Island populations share a more recent history than others, but were subsequently isolated from one another during climatic warming. Together, these results show that M. emersoni populations have responded directly to climate warming by adapting physiologically and modifying their distribution.

Fig. 2
figure 2

M. emersoni populations have adapted physiologically to the prevalent environmental conditions that were established after climate warming began 10,000 years BP. a Heat shock experiment: M. emersoni workers from warmer collection sites show significantly less mortality after a one hour heatshock at 50 ° C than workers from colder collection sites (R2 = 0.42, p = 0.035), indicating that M. emersoni colonies are adapted to the local environmental conditions that were established after climatic warming. b Cold shock experiment: M. emersoni workers from colder collection sites show a significantly more rapid recovery after a 16h cold shock at 4 ° C than workers from warmer collection sites (R2 = 0.46, p = 0.018)

Fig. 3
figure 3

M. emersoni populations have been isolated from each other over the last 200,000 years and modern population structure reflects geography. a Simulations showing approximate timing of isolation events between M. emersoni populations in the Arizona Sky Islands during the last 200,000 years (see also Additional file 1: Figure S2). These simulations revealed two geographical groupings of Sky Islands: North and South. b Redundancy analysis showing that most of the mtDNA genetic variation is explained by geography (site). The Venn diagram represents the percentage of the genetic variation that can be attributed to each factor (circles) or their combined effect (overlapping areas). This analysis supports the hypothesis of multiple origins of the wingless phenotype, where genetic variation is structured by geography rather than habitat, phenotype or elevation (see also Additional file 1: Figure S3 and S4)

We discovered two alternative queen phenotypes, winged and wingless, which likely reflect different life histories (independent and dependent colony foundation) adapted to different habitats. We therefore determined whether the frequency of winged and wingless phenotypes changes along the ecological gradients on each Sky Island. We found that wingless queens do indeed differ in their frequencies along the ecological gradient (Additional file 1: Figure S4 and S5), suggesting that there may be a set of environmental factors that the winged and wingless phenotypes may be adapting to on each Sky Island. To potentially explain the distribution of winged and wingless queens across collection sites, we used previous literature to construct 12 biological models (Table 1). Each model included one or several of 5 key environmental factors as well as possible interactions between those factors (Additional file 1: Table S1). We tested the ability of the 12 competing models to describe the data using a model-ranking approach [51]. Model m06, which includes the combined effects of habitat fragmentation and elevation, is by far the model with the highest probability (0.58) accounting for 50 % of the variability in the data (Table 1, Additional file 1: Table S2). Four out of the 5 highest ranking models included habitat fragmentation and/or elevation as explanatory variables, bringing additional support for a combined effect of these variables on the distribution of wingless queens along the ecological clines. Indeed, significant habitat fragmentation has been postulated to favor the evolution of wingless queens for short distance dispersal (dependent colony foundation) over winged queens for long distance dispersal (independent colony foundation), because winged queens have increased chances of landing in an inhospitable habitat. These results support a scenario where climate warming indirectly drove the evolution of the wingless queen phenotype through habitat changes. The paleoecological record indicates that during the last glaciation 90,000 to 20,000 years BP forests were ubiquitous and continuously distributed at the lower elevations [24, 25]. Because winged queens occur more frequently in low elevation and more continuous habitats, we can infer that the ancestral populations were likely composed of winged queens. This is consistent with our finding of significant gene flow between all five Sky Islands 90,000 years BP (Additional file 1: Figure S1). Subsequently, wingless queens emerged when climatic warming drove the appearance of fragmented high elevation habitats on each of the Sky Islands 20,000 years BP.

Table 1 A combination of landscape fragmentation and elevation best explain the distribution of wingless queens along ecological gradients. Based on previous studies, we generated 12 models (M01 to M12) representing alternative hypotheses to explain the variation in frequency of wingless queens along the ecological gradient of each Sky Islands. The 12 models were ranked using the method of [121] implemented in R (Additional file 1: Table S2) and model probabilities (Prob) are reported, with m06 scoring as the highest ranking model

Our data show that emergence of the wingless queen phenotype occurred repeatedly by parallel evolution in response to similar ecological conditions arising on each of the Sky Islands during climatic warming. First, modern gene flow between mountain ranges is absent or low; the degree of genetic differentiation (F ST ) is high between collection sites on different Sky Islands, but is lower within a Sky Island (Additional file 1: Table S3 and S5). Second, phylogenetic and redundancy analyses based respectively on multilocus AFLP data and mtDNA sequences show that M. emersoni colonized each Sky Island independently (Fig. 3 b, Additional file 1: Figure S3 and S4). We find that populations on each Sky Island contain unique haplotype variants, while no shared haplotypes were recovered among Sky Islands, indicating that unique haplotypes originated within each Sky Island (Additional file 1: Figure S3). Finally, our results show that geography dictates the patterns of genetic differentiation at neutral loci: winged and wingless queens from the same Sky Island share a more recent ancestor than queens from different Sky Islands (Additional file 1: Figure S3 and S4), and genetic variation among queens is best explained by geography (collection site) rather than habitat or phenotype (Fig. 3 b). All together, this supports an incipient process of parallel differentiation of the populations and suggests that M. emersoni wingless queens emerged independently on each Sky Island when each population encountered similar but newly established ecological conditions, namely, fragmented and high elevation habitats.

In many taxa, the repeated evolution of a trait observed across species within a genus, such as winglessness in stick insects [52], Drosophila wing spots [53], or armor plate loss in freshwater sticklebacks [1], have not only been generated by similar selective pressures but also by ancestral developmental potentials [3, 5456]. Developmental potentials are unexpressed genetic programs, such as gene networks, or unexpressed developmental thresholds that determine sensitivity of these networks to specific environmental conditions. Developmental potentials are ancestral because they are often retained for millions of years in modern groups of organisms, most probably because of pleiotropy [3, 57]. In ants, Rajakumar et al. [54] demonstrated that ancestral developmental potentials can be induced by mutations or environmental perturbations to produce ancestral phenotypic variants (anomalies) that selection acts on to facilitate parallel evolution. We therefore searched the literature for evidence of an ancestral developmental potential to produce wingless queens in the genus Monomorium. We found that, unlike other genera, a large proportion of New World species (10 out of 16) in Monomorium show variable frequencies of appearance of wingless queens (Table 2). In 5 species (M. emersoni, M. ergatogyna, M. trageri, M. viridae and in M. sp. AZ-03), queens are winged but the existence of wingless phenotypes have been reported either as anomalies (sporadic occurrences) or alternative phenotypes (common occurrence of both phenotypes), while in the other five species (M. compressum, M. cyaneum, M. ebeninum, M. sp. AZ-01 and M. sp. cf. ergatogyna), wingless queens appear as a fixed phenotype in the colony. The 6 remaining species are known to have winged queens only. The common occurrence of wingless queens as sporadic anomalies in colonies, as well as the evolution of species with alternative or fixed phenotypes, indicates the existence of an ancestral developmental potential to produce wingless queens in Monomorium, which likely played a role in facilitating the parallel evolution of the wingless queen phenotype within M. emersoni. This developmental potential to produce wingless queens could be unique to the genus Monomorium, or may be a more ancient feature that facilitated the 200 independent evolutionary origins of wingless queens across ants. The ancestral developmental potential to produce wingless queens may have been co-opted from the wing polyphenism in ants - the ability of a single genome to produce winged queens and wingless workers in a colony in response to environmental cues [58], just after the origin of ants. It is unlikely, however, that wingless queens evolved by directly co-opting this capacity to produce the wingless worker caste de novo during the evolution of M. emersoni, because wingless workers in M. emersoni interrupt wing development much earlier in development than wingless queens (Additional file 1: Figure S5). Collectively, our results suggest that in response to climatic change, selection on the phenotypic expression of the ancestral potential for producing wingless queens may have driven the parallel evolution of the wingless phenotype on each Sky Island.

Table 2 The presence of wing loss in queens in Monomorium minimum-group species. Summary of known queen phenotypes for Monomorium minimum-group species that are native to the United States and Mexico; species are listed alphabetically

Based on this shared ancestral developmental potential for producing wingless queens in the genus Monomorium, as well as the very recent history shared by M. emersoni populations, we initially predicted that similar changes in the wing developmental system would underlie the parallel evolution of the wingless phenotype. Wings in holometabolous insects develop from imaginal discs that are clusters of cells in larvae [17]. These cells also give rise to the flight muscles and to the wing hinge during metamorphosis [20, 5962]. In wingless queens, wing imaginal discs appear as vestigial remnants during the last larval instar (Additional file 1: Figure S5). We characterized the cell division and growth of wing imaginal discs and found that vestigial wing discs of wingless queen larvae from different mountain ranges significantly differ in size (Fig. 4) and levels of cell division (Fig. 4 b to f). Because larvae were reared under the same experimental conditions, these differences are a consequence of an evolutionary change in growth, not developmental plasticity, in wingless queens across Sky Islands. Surprisingly, we discovered that these differences may be the consequence of a demographic split that happened 80,000 years BP between the northern (Catalinas and the Pinals) and the southern (Huachucas and the Chiricahuas) Sky Islands (Figs. 1 e, 3 a and 4 g). Wingless queens from northern Sky Islands possess smaller discs and low levels of cell division (Fig. 4 h), whereas those from the southern Sky Islands possess large imaginal discs and high levels of cell division. This association suggests that demographic events drove the evolution of novel variations in wing disc growth across M. emersoni populations.

Fig. 4
figure 4

Wingless queens from the Northern and the Southern Sky Islands show different growth patterns. a Surface area (μm2) of leg discs (x-axis) versus vestigial wing disc (y-axis) for wingless queen last instar larvae from different mountain ranges. A line was fitted through with a Standardized Major Axis bivariate line-fitting analysis. The slopes are significantly different from each other (p-value = 0.00017), reflecting an evolutionary change affecting developmental trajectories of imaginal discs among different mountains since their isolation. b Winged queen larvae shows extensive cell division in imaginal discs in last larval instar. Wingless queen larvae from the Northern Sky Islands c Pinals and the d Catalinas show very few cell divisions in their vestigial wing discs while the Southern Sky Islands, the e Huachucas and f Chiricahuas show more. Red is PH3, green is DAPI, and yellow arrowheads indicate a cell that is expressing PH3. g Map of the Sky Islands indicating the position of Northern and Southern Sky Islands. h Graph of the number of dividing cells per unit area in vestigial discs of wingless queen larva from different Sky Islands. The x-axis shows the different mountain ranges and the y-axis the number of cells dividing per (μm2). The number of dividing cells per unit area in vestigial wing discs differs between Northern and Southern Sky Islands (ANOVA, p <0.0001). Thus, the vestigial wing discs of queens from Southern Sky Islands show significantly higher levels of cell division and are larger than vestigial wing discs of queens from Northern Sky Islands

Wing discs express a gene network that controls growth and patterning of wings (Fig. 1 f), of wing muscles, and of the wing hinge (Fig. 7, last column), and has been highly conserved across holometabolous insects for approximately 325 million years [20, 6366]. In ants, the expression of this network is conserved in winged queens and males but is interrupted at specific points in the network in wingless workers [63]. These points of interruption in wingless workers are spatiotemporally associated with patterns of apoptosis in the vestigial wing disc [65]. We therefore determined whether points of interruption in this network in wingless queens has evolved in parallel or whether it differs between Sky Island populations. To confirm that expression of genes in this network is conserved in M. emersoni winged queens, we selected 6 genes within this network (Fig. 1 f), including 4 genes involved in patterning wing tissue (engrailed [en], Ultrabithorax [ubx], cut [cut], wingless [wg]), one gene involved in wing muscle precursor development (myosin enhancer factor 2 [mef2]), and one gene involved in wing hinge development (extradenticle [exd]). In winged queens, the expression of these genes is similar to those in the wing discs of Drosophila [20, 5961] and winged queens and males of other ant species [63, 67] (Fig. 5 b, c, and d, Fig. 6 b, Fig. 7 b and Additional file 1: Figure S6B), indicating that this network is conserved.

Fig. 5
figure 5

Expression profiles of cut, Ubx, and en in vestigial discs of wingless queens from different Sky Islands show similar expression patterns. Panels in the left column a, e, i, m, and q represent diagrams of expression patterns of Cut (blue), Ubx (violet) and En (pink). From left to right, columns show gene expression profiles (purple) of Cut (b, f, j, n, r), Ubx c, g, k, o, s, and En (d, h, l, p, t) in winged queens and wingless queens across four Sky Islands. Black arrowheads point to vestigial discs and asterisks indicate an absence of expression in the vestigial wing discs. Insets show forewing discs at higher magnification

Fig. 6
figure 6

Wingless queen larvae develop wing muscle precursors, but lack flight muscles entirely. Panels on the left column (a, e, i, m, q) represent diagrams of Mef-2 expression. Mef-2 expression in larval myoblasts is shown as a localized stain in nuclei of myoblasts associated with the imaginal discs b, f, j, n, and r. Details on Mef-2 expression are given in Additional file 1: Figure S7. Inserts show forewing discs at higher magnification. Dashed lines on Mef2 expression patterns highlight the vestigial disc boundaries based on DAPI staining. Virtual sections (microCT) through the thorax (last two columns) showing the internal anatomy of a winged queen (top row) and wingless queens from all Sky Islands. c, g, k, o, and s show transverse sections and d, h, l, p, and t show sagittal sections. Longitudinal muscles are artificially colored in pink and dorso-ventral muscles in blue. Wingless queens from all Sky Islands show a complete absence of indirect flight muscles. These results show that, regardless of their origin, wingless queens do not develop the main flight muscles

Fig. 7
figure 7

The variation in Exd expression pattern across Sky Islands corresponds to changes of fine details in thorax morphology of wingless queens. Panels on the left column (a, d, g, j, m) represent diagrams of Exd expression. Exd expression in winged queens b and wingless queens e, h, k, and n. Black arrowheads point to Exd expression in imaginal discs. Insets show a close-up of Exd expression in the vestigial imaginal discs. 3D X-ray microtomography thoracic morphology of a winged queen c and wingless queens from Pinals f, Catalinas i, Huachucas l, and Chiricahuas o Sky Islands. Insets show a close-up of the wing hinge region. Yellow brackets indicate vestigial forewing hinge and yellow arrowheads indicate outgrowths of the hindwing hinge. Note that the vestigial forewing hinge in f and o is prominent with two distinct outgrowths, while it is reduced in i, l with only one prominent outgrowth. The vestigial hindwing hinge is present in all populations

We then determined whether the expression of these genes has been interrupted in wingless queens, and found that their expression is altered relative to winged queens (Figs. 5, 6 and 7 and Additional file 1: Figure S6). In wingless queens from different Sky Islands, we found that the 4 genes involved in wing tissue patterning (ubx, en, cut, wg) show similar aberrant patterns of expression (Fig. 5 e to t and Additional file 1: Figure S6C to H). These shared alterations reveal that wingless queens, regardless of their geographic origin and extent of vestigial wing disc growth (Fig. 4 a), interrupt wing patterning at least in part through similar changes in gene expression. Therefore, although organ growth retains the signature of historical events caused by past demographic changes, alteration of the underlying expression of wing patterning genes is recurrent among populations on each Sky Island.

In contrast to genes involved in wing patterning, we observed that the expression of mef2 - a key transcription factor involved in wing muscle development [61] - differs among populations and shows a signature of the demographic split between the southern and the northern Sky Islands. Mef2, which is expressed in the nuclei of myoblasts and appears as a dotted pattern in wing discs, shows expression differences between wingless queens from northern and southern Sky Islands (Fig. 6 f and j versus N and R, and Additional file 1: Figure S7). In the northern Sky Islands, Mef2 expression is restricted to the outer edge of the vestigial disc in wingless queens, whereas in the southern Sky Islands, it is expressed more extensively and covers the central region of vestigial discs. Thus, like patterns of wing disc growth (Fig. 4 h), Mef2 expression patterns precisely recapitulate the history of contacts between Sky Islands. Furthermore, these differences in the pattern of Mef2 expression in vestigial wing discs are likely to have evolved neutrally because wingless queens lack flight muscles regardless of their origin (Fig. 6 g, h, k, l, o, p, s and t). This suggests that another mechanism, such as the interruption of other genes in the network or apoptosis of myoblasts [65], prevents mef2 expressing cells from further differentiating into flight muscles.

Finally, we found that exd - a gene known to be responsible for wing hinge development [20] - shows highly distinct expression patterns among wingless queens from different Sky Islands (Fig. 7 b, e, h, k, and n) that correspond to fine differences in wing hinge morphology (Fig. 7 c, f, i, l, and o). In wingless queens coming from the same Sky Island, the aberrant pattern of Exd expression is the same across individuals (Additional file 1: Figure S8). However, in wingless queens from the Pinals (north) and the Chiricahuas (south), Exd is expressed over the entire vestigial disc and adults have a prominent vestigial fore- and hindwing hinges, whereas in wingless queens from the Catalinas (north) and Huachucas (south), Exd is only expressed in the anterior edge of vestigial discs and adults have a prominent vestigial hindwing hinge but the forewing hinge is absent or dramatically reduced in size (Fig. 7 and Additional file 1: Table S6). These observed differences within northern (Pinals and Catalinas) and within southern (Chiricahuas and Huachucas) Sky Islands indicate that these differences arose following the complete isolation of the five Sky Island populations 10,000 years BP and are therefore an indirect consequence of climatic warming (Fig. 8).

Fig. 8
figure 8

Summary of climatic changes, population history, wingless phenotype evolution and the associated changes in developmental systems in M. emersoni. The population tree represents the history of contacts during the last 200,000 years. Ancestral potential for wing loss (grey ant profile). Climatic changes are depicted in the bar on the left, with the associated changes in the habitats in different colors (green: largely continuous landscape, yellow: fragmented landscape). Changes in the vestigial disc growth, wing hinge morphology, wing muscles and in expression patterns of genes responsible for wing, wing hinge and wing muscle development are depicted by symbols on the tree. Gene expression on vestigial discs is depicted by purple (Exd), dotted green (Mef2) or black (wing patterning genes: Ubx, En, Cut and wg). Changes in cell division are depicted by the red dots


A paleoecological reconstruction for the parallel evolution of a developmental network and its phenotype in response to past climate change on Sky Islands

We have provided evidence for the following scenario, which shows that in response to both direct and indirect effects of past climate change on Sky Islands, a developmental gene network and its corresponding phenotype can respond to changing ecological conditions on short time scales (Fig. 8): Paleoecological reconstructions of the Arizona Sky Islands [24, 25] show that during the last glaciation between about 90,000 to 20,000 years BP, a continuous forest habitat across the landscape prevailed, except for high elevations which were inhospitable alpine habitats. The developmental potential to produce wingless queens in M. emersoni was present during this period. During or immediately following the vicariant event 80,000 years BP causing the demographic split between northern and southern Sky Islands, these M. emersoni populations accumulated unexpressed differences in the genes that control wing growth and flight muscle development. Subsequently, during climatic warming 10,000 years BP, the landscape became fragmented and new habitats arose at high elevations that were suitable for colonization by M. emersoni. As populations encountered these new habitats on each of the Sky Islands, the wingless phenotype evolved independently and became more frequent through selection on the phenotypic expression of the developmental potential for producing wingless queens. The parallel evolution of wingless queens on each Sky Island occurred through recurrent interruptions of similar wing-patterning genes in the network (ubx, en, cut, wg), and at the same time, induced the release of the unexpressed differences in vestigial wing growth and mef-2 expression that accumulated in populations during or after the demographic split between the northern and southern Sky Islands 80,000 years BP. Climate warming during the last 10,000 years BP formed vast deserts between Sky Islands thereby isolating the five Sky Island populations from one another and forced M. emersoni populations to adapt physiologically to the warm conditions at the lower elevations of their range. This may have led to the rapid evolution of the vestigial wing hinge and its underlying exd expression on each Sky Island.

Novelty as a mosaic of contingency and determinism

We show that developmental systems can facilitate two types of evolutionary change that contribute to the differentiation of populations within a species: one that is governed by chance events, is unpredictable, and has produced small differences in phenotype and genotype: wing hinge morphology, wing disc growth and expression of exd and mef-2. In contrast, the other type is repeatable and has produced large differences in phenotype and genotype: the wingless queen. The evolution of wingless queens represents an evolutionary novelty - it is the gain of a developmental switch, the complete loss of wings and wing muscles, multiple alterations in the expression of the core developmental network underlying wing development, and the gain of a new life history strategy (dependent colony foundation). These findings offer a new perspective for the ongoing discussion of novelty generation in phenotypic evolution [68]. On the one hand, it has been argued that adaptive solutions through natural selection are deterministic and that evolution is largely predictable [69], while on the other, it has been argued that historical contingency, which is the unique chain of prior events that determine which of several evolution paths a lineage will follow, renders evolution unpredictable [70]. Several studies, from a phylogenetically diverse array of organisms, show that the parallel evolution of phenotypic traits often occurs through recurrent alterations in the same genes [12, 23, 7174]. Yet, at the same time, these and other studies report slight differences across multiple levels of biological organization: phenotype, development, and genetics [7577]. In vitro studies on microbial populations suggest that these similarities and differences underlying the parallel evolution of traits are the consequence of a complex interplay between deterministic processes and historical events [7880]. We show, in natural populations within a species, that the emergence of novelty at multiple biological levels can be generated by the combination of ancestral developmental potentials facilitating parallel evolution in response to similar ecological conditions and chance events that create slight variations on these repeatable evolutionary novelties.

A mosaic of contingent and deterministic processes at multiple biological levels may be a general property of the evolution of novelty. Most, if not all, organisms experience similar population dynamics over evolutionary time and parallel evolution within and between species is a nearly universal pattern throughout the tree of life [69, 81]. Many independently evolved traits, such as, pigmentation spots on Drosophila wings [53, 82], pharyngeal jaws of cichlid fish [83], hind limbs of anoles lizards [84], coat color of wild mice [85], and coloration on butterfly wings [86], exhibit patterns consistent with our findings in that they show slight variations in phenotype and genotype. Our results suggest that these slight variations are generated by microevolutionary forces and historical contingency acting between populations, together causing the gradual accumulation and evolution of unexpressed changes in regulatory pathways. These unexpressed differences become expressed during the process of parallel evolution leading to novel variations on independently evolved traits. This process driving parallel phenotypic novelty within species may translate into parallel phenotypic novelty between species, and as in the case pharyngeal jaws of cichlid fish [83], slight variations accumulated in one of the lineages may potentially fuel phenotypic diversification.

Implications for bridging Micro- and Macroevolution

One of the biggest challenges facing evolutionary theory today is to explain how microevolution translates into macroevolution. Evolutionary biologists have addressed this problem from the perspective of paleontology [6, 87], population and quantitative genetics [88], and evolutionary developmental biology [16, 89, 90], producing a range of different perspectives on the translation between micro- and macroevolution. One perspective, which took hold during the neo-Darwinian synthesis in the 1940’s, proposes no distinction micro- and macroevolution, where macroevolution is simply the product of a specific kind of microevolution (mutations of small effect, individual level selection and drift) extrapolated over long temporal scales [88]. This perspective, however, is not well supported in the face of recent advances in quantitative genetics [8, 91], experimental evolution [92], developmental plasticity [3] and paleontology [93]. These recent advances show that: (1) dynamics of evolutionary change observed between individuals within populations cannot be extrapolated to explain the dynamics of change observed between species within clades [87, 93]; and (2) alleles of large effect [8, 91] and large discrete phenotypic variants [3, 54, 94] can be fixed in populations. Indeed, most modern textbooks of evolution now distinguish micro- from macroevolution and use the species boundary as the dividing line between them, most often defining microevolution as ‘evolution below the species level’ - small changes, adaptive or neutral, that arise within populations, while macroevolution is most often defined as ‘evolution above the species level’ - large changes that arise within populations leading to the origin of species or higher taxonomic groups [2, 9598]. Although all changes must necessarily arise within populations, macroevolutionary changes are thought to arise during or after populations have undergone speciation, and therefore, characterize differences observed between closely-related species or higher taxonomic groups [3, 93].

Our findings, however, suggest that the species boundary may not adequately separate micro- and macroevolutionary changes. We observed macroevolution-like change (emergence of wingless queens) between populations within the same species (M. emersoni), which is comparable to the variation observed between species in the genus Monomorium (Table 2). Our observations are consistent with an idea known as ‘intraspecific macroevolution’, which proposes that micro- and macroevolution are distinct, but that speciation is not the primary cause of large phenotypic changes between species leading to macroevolution [3]. It proposes that both small and large phenotypic changes occur within the same species, and that the existence of alternative phenotypes within lineages facilitates macroevolution by reducing the negative fitness consequences associated with the appearance and fixation of large phenotypic variants [3].

Finally, it has been proposed that changes in hierarchically organized developmental networks can simultaneously produce the large- and small-scale changes associated with macro- and microevolution [99, 100]. Genes in the network involved in core developmental processes (‘kernels’ sensu Davidson, [18]), such as body plan formation, are thought to produce large-scale variants when changed, be highly stable during evolution, and are associated with macroevolution. In contrast, genes in the network involved in terminal differentiation processes, like pigment formation, are expected to produce small-scale variants, be more labile during evolution, and are associated with microevolution [99, 100]. Our results are consistent but show that the core developmental network responsible for wing formation facilitates the evolution of both types of changes within a species in response to ecological change. Based on all of our findings, we propose that the parallel phenotypic novelty within species may translate into parallel phenotypic novelty between species providing a possible bridge between micro- and macroevolution.

The role of organismal development in past and future climate change

Understanding how developmental systems respond to past climate change is critical to facing the challenges of global warming over the next few decades. Integrating the responses of developmental systems with environmental changes will be necessary for building accurate projection models to predict species distributions under climatic change scenarios. Predicting species distributions accurately can affect decision making about conservation biology strategies, agriculture and public health issues. Current projection models integrate demographic processes such as migration and population size changes, but few integrate the possible response of developmental systems to environmental changes, which as we show, may profoundly affect the life history and dispersion abilities of a species. Ultimately, such changes directly or indirectly affect species interactions and ecosystem processes [101]. Sky Island mountain forests have been used as models and indicators of climate change. They have responded to past [24] and present [36] global warming by contracting towards the higher elevations [35], and climatic and ecological models predict that warming and desertification in Arizona will continue over the next 50 years [33, 34, 102]. Several species found in the mountainous areas of the American Southwest already show signs of changes in phenology in response to recent global warming [103, 104]. These dramatic changes in environmental conditions will be integrated by the developmental systems of different organisms, and as we show, can generate both predictable and unpredictable phenotypic responses. Generally integrating the predictable part of the response of developmental systems into climate change studies should hold the promise of making more accurate predictions for organismal responses to climate change.


Since the modern synthesis, our view of how population differentiation occurs in the wild is largely based on the advances made in population and quantitative genetics, especially with the advent and integration of genomics tools. This view is one of changes in allele and phenotypic frequencies in response to adaptive and non-adaptive forces. By integrating organismal development and paleoecology on Sky Islands, we enrich this view by showing that the emergence of novelty within populations occurs through a combination contingent and deterministic processes. This combination results in the parallel evolution of traits with slight variations at multiple levels of biological organization (phenotype, growth, developmental networks). We propose that this parallel phenotypic novelty within the species may represent the springboards for translation into macroevolutionary changes between species. Studying these evolutionary processes on the Sky Islands, with their rich fauna provides replicated ecological gradients, are a good model system for understanding how development influences an organism’s response to past and future climate change.

Materials and methods

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. 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.

Data availability

The population sets of DNA sequences generated in this study are available on GenBank (accession numbers KT356282-KT356545). Other data used in this study are available upon request to the authors.


  1. Schluter D, Nagel L. Parallel speciation by natural selection. American Naturalist. 1995; 146:292–301.

    Article  Google Scholar 

  2. Futuyma DJ. Evolutionary Biology: Sinauer; 1998.

  3. Eberhard MJW. Developmental Plasticity and Evolution. New-York: Oxford University Press; 2003.

    Google Scholar 

  4. Darwin C. On the Origins of Species by Means of Natural Selection. London: Murray; 1859.

    Google Scholar 

  5. Mayr E. Systematics and the Origin of Species, from the Viewpoint of a Zoologist: Harvard University Press; 1942.

  6. Gould SJ, Eldredge N. Punctuated equilibria: an alternative to phyletic gradualism In: Schopf TJM, editor. Models in Paleobiology. San Francisco: Freeman Cooper: 1972. p. 82–115.

    Google Scholar 

  7. Raff RA. The Shape of Life: Genes, Development, and the Evolution of Animal Form. Chicago: University of Chicago Press; 1996.

    Google Scholar 

  8. Orr HA. The population genetics of adaptation: the distribution of factors fixed during adaptive evolution. Evolution. 1998;:935–49.

  9. Hartl DL, Clark AG, Clark AG. Principles of Population Genetics, vol. 116. Sunderland: Sinauer associates; 1997.

    Google Scholar 

  10. Nielsen R. Molecular signatures of natural selection. Annu. Rev. Genet. 2005; 39:197–218.

    Article  CAS  PubMed  Google Scholar 

  11. Jones FC, Chan YF, Schmutz J, Grimwood J, Brady SD, Southwick AM, et al. A genome-wide snp genotyping array reveals patterns of global and repeated species-pair divergence in sticklebacks. Current Biology. 2011; 22(1).

  12. Chan Y, Marks M, Jones F, Villarreal G, Shapiro M, Brady S, et al. Adaptive evolution of pelvic reduction in sticklebacks by recurrent deletion of a pitx1 enhancer. Science. 2010; 327.5963:302–5.

    Article  CAS  Google Scholar 

  13. Arnegard ME, McGee MD, Matthews B, Marchinko KB, Conte GL, Kabir S, et al. Genetics of ecological divergence during speciation. Nature. 2014.

  14. Soria-Carrasco V, Gompert Z, Comeault AA, Farkas TE, Parchman TL, Johnston JS, et al. Stick insect genomes reveal natural selection’s role in parallel speciation. Science. 2014; 344(6185):738–42.

    Article  CAS  PubMed  Google Scholar 

  15. Pfennig DW, Wund MA, Snell-Rood EC, Cruickshank T, Schlichting CD, Moczek AP. Phenotypic plasticity’s impacts on diversification and speciation. Trends in Ecology & Evolution. 2010; 25(8):459–67.

    Article  Google Scholar 

  16. Nunes MD, Arif S, Schlötterer C, McGregor AP. A perspective on micro-evo-devo: progress and potential. Genetics. 2013; 195(3):625–34.

    Article  PubMed Central  PubMed  Google Scholar 

  17. Gilbert SF. Developmental Biology. Sunderland: Sinauer Associates; 2010.

    Google Scholar 

  18. Davidson EH. The Regulatory Genome: Gene Regulatory Networks in Development and Evolution. Burlington: Academic Press; 2010.

    Google Scholar 

  19. Gilbert SF, Epel D. Ecological Developmental Biology. Sunderland: Sinauer Associates; 2009.

    Google Scholar 

  20. Carroll SB, Grenier JK, Weatherbee SD. From DNA to Diversity: Molecular Genetics and the Evolution of Animal Design. Cambridge: Blackwell; 2001.

    Google Scholar 

  21. Marigo V, Scott MP, Johnson RL, Goodrich LV, Tabin CJ. Conservation in hedgehog signaling: induction of a chicken patched homolog by sonic hedgehog in the developing limb. Development. 1996; 122(4):1225–33.

    CAS  PubMed  Google Scholar 

  22. Stern DL. Evolution, Development, & the Predictable Genome: Roberts and Company Publishers; 2011.

  23. Stern D, Orgogozo V. Is genetic evolution predictable?Science. 2009; 323(5915):746–51.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Lomolino M, Brown J, Davis R. Islands biogeography of montane forest mammals in the american southwest. Ecology. 1989;:180–94.

  25. Thompson RS, Anderson KH. Biomes of western North America at 18,000, 6000 and 0 C-14 yr bp reconstructed from pollen and packrat midden data. Journal of Biogeography. 2000; 27(3):555–84.

    Article  Google Scholar 

  26. Warshall P. The madrean sky island archipelago: a planetary overview In: DeBano L, Ffolliott P, Ortega-Rubio A, Gottfried G, Hamre R, Edminster C, editors. RM-GTR-264 : Biodiversity and Management of the Madrean Archipelago: the Sky Islands of Southwestern United States and Northwestern Mexico. Fort Collins: US Department of Agriculture, Forest Service: 1994. p. 6–18.

    Google Scholar 

  27. Gottfried GJ, Gebow BS, Eskew LG, Edminster CB. Connecting Mountain Islands and Desert Seas: Biodiversity and Management of the Madrean Archipelago II. Proc. RMRS-P-36. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station: 1–6.

  28. Ober K, Matthews B, Ferrieri A, Kuhn S. The evolution and age of populations of Scaphinotus petersi Roeschke on Arizona Sky Islands (coleoptera, carabidae, cychrini). ZooKeys. 2011; 147:183.

    Article  PubMed  Google Scholar 

  29. Heinze J, Rueppell O. The frequency of multi-queen colonies increases with altitude in a nearctic ant. Ecological Entomology. 2014;39:4.

  30. Masta S. Phylogeography of the jumping spider Habronattus pugillis (Araneae : Salticidae): Recent vicariance of sky island populations?Evolution. 2000; 54(5):1699–711.

    Article  CAS  PubMed  Google Scholar 

  31. Smith CI, Farrell BD. Range expansions in the flightless longhorn cactus beetles, Moneilema gigas and Moneilema armatum, in response to pleistocene climate changes. Molecular Ecology. 2005; 14(4):1025–44.

    Article  CAS  PubMed  Google Scholar 

  32. Perez-Alquicira J, Molina-Freaner F, Pinero D, Weller S, Martinez-Meyer E, Rozas J, et al. The role of historical factors and natural selection in the evolution of breeding systems of Oxalis alpina in the sonoran desert sky islands. Journal of Evolutionary Biology. 2010; 23(10):2163–175.

    Article  CAS  PubMed  Google Scholar 

  33. Kupfer JA, John A, Cairns DM. The suitability of montane ecotones as indicators of global climatic change. Progress in Physical Geography. 1996; 20(3):253–72.

    Article  Google Scholar 

  34. Kupfer JA, Balmat J, Smith JL. Shifts in the potential distribution of sky island plant communities in response to climate change. USDA Forest Service Proceedings. 2005;:485–90.

  35. Brusca RC, Wiens JF, Meyer WM, Eble J, Franklin K, Overpeck JT, et al. Dramatic response to climate change in the southwest: Robert Whittaker’s 1963 Arizona mountain plant transect revisited. Ecology and Evolution. 2013; 3(10):3307–319.

    PubMed Central  PubMed  Google Scholar 

  36. Garfin G, Jardine A, Merideth R, Black M, LeRoy S. Assessment of Climate Change in the Southwest United States: Island Press; 2013.

  37. Hölldobler B, Wilson EO. The Ants. Cambridge: The Belknap Press of Harvard University Press; 1990.

    Book  Google Scholar 

  38. Andersen AN, Majer JD. Ants show the way down under: invertebrates as bioindicators in land management. Frontiers in Ecology and the Environment. 2004; 2(6):291–8.

    Article  Google Scholar 

  39. Agosti D. Ants: Standard Methods for Measuring and Monitoring Biodiversity. Biological diversity handbook series: Smithsonian Institution Press; 2000.

  40. Diamond SE, Nichols LM, McCoy N, Hirsch C, Pelini SL, Sanders NJ, et al. A physiological trait-based approach to predicting the responses of species to experimental climate warming. Ecology. 2012; 93(11):2305–312.

    Article  PubMed  Google Scholar 

  41. Pelini SL, Diamond SE, MacLean H, Ellison AM, Gotelli NJ, Sanders NJ, et al. Common garden experiments reveal uncommon responses across temperatures, locations, and species of ants. Ecology and Evolution. 2012; 2(12):3009–015.

    Article  PubMed Central  PubMed  Google Scholar 

  42. Hölldobler B, Wilson EO. The number of queens: an important trait in ant evolution. Naturwissenschaften. 1977; 64(1):8–15.

    Article  Google Scholar 

  43. Heinze J, Keller L. Alternative reproductive strategies: a queen perspective in ants. Trends in Ecology & Evolution. 2000; 15(12):508–12.

    Article  Google Scholar 

  44. Peeters C. Convergent evolution of wingless queens reproductive across all subfamilies of ants, and sporadic loss of winged queens (hymenoptera: Formicidae). Myrmecological News. 2012; 16:75–91.

    Google Scholar 

  45. Peeters C, Ito F. Colony dispersal and the evolution of queen morphology in social hymenoptera. Annual Review of Entomology. 2001; 46(1):601–30.

    Article  CAS  PubMed  Google Scholar 

  46. Heinze J, Buschinger A. Queen polymorphism in Leptothorax Species-A - its genetic and ecological background (hymenoptera, formicidae). Insectes Sociaux. 1989; 36(2):139–55.

    Article  Google Scholar 

  47. Roff DA. The evolution of flightlessness in insects. Ecological Monographs. 1990;:389–421.

  48. Favé M-J. The role of past climate change in driving novelties in sky island ant populations. Montreal: PhD thesis, McGill University; 2013.

    Google Scholar 

  49. Cremer S, Heinze J. Stress grows wings: Environmental induction of winged dispersal males in cardiocondyla ants. Current Biology. 2003; 13(3):219–223.

    Article  CAS  PubMed  Google Scholar 

  50. Mitchell HK, Petersen NS. Developmental abnormalities in drosophila induced by heat shock. Developmental genetics. 1982; 3(2):91–102.

    Article  CAS  Google Scholar 

  51. Whittingham MJ, Stephens PA, Bradbury RB, Freckleton RP. Why do we still use stepwise modelling in ecology and behaviour?Journal of animal ecology. 2006; 75(5):1182–9.

    Article  PubMed  Google Scholar 

  52. Whiting MF, Bradler S, Maxwell T. Loss and recovery of wings in stick insects. Nature. 2003; 421(6920):264–7.

    Article  CAS  PubMed  Google Scholar 

  53. Prud’Homme B, Gompel N, Rokas A, Kassner VA, Williams TM, Yeh SD, et al. Repeated morphological evolution through cis-regulatory changes in a pleiotropic gene. Nature. 2006; 440(7087):1050–3.

    Article  PubMed  CAS  Google Scholar 

  54. Rajakumar R, Mauro D, Dijkstra M, Huang M, Wheeler D, Hiou-Tim F et al. Ancestral developmental potential facilitates parallel evolution in ants. Science. 2012; 335(6064):79–82.

    Article  CAS  PubMed  Google Scholar 

  55. Wund MA, Valena S, Wood S, Baker JA. Ancestral plasticity and allometry in threespine stickleback reveal phenotypes associated with derived, freshwater ecotypes. Biological Journal of the Linnean Society. 2012; 105(3):573–83.

    Article  PubMed Central  PubMed  Google Scholar 

  56. Losos JB. Convergence, adaptation, and constraint. Evolution. 2011; 65(7):1827–40.

    Article  PubMed  Google Scholar 

  57. Abouheif E, Favé M-J, Ibarrarán-Viniegra AS, Lesoway MP, Rafiqi AM, Rajakumar R. Eco-evo-devo: The time has come. In: Ecological Genomics. Springer: 2014. p. 107–25.

  58. Molet M, Wheeler DE, Peeters C. Evolution of novel mosaic castes in ants: Modularity, phenotypic plasticity, and colonial buffering. The American Naturalist. 2012; 180(3):328–41.

    Article  PubMed  Google Scholar 

  59. (Cohen SM, Bate M, Arias AM, editors.)1993. The Development of Drosophila Melanogaster vol. 1. Cold Spring Harbor Laboratory Pr.

  60. Weatherbee SD, Halder G, Kim J, Hudson A, Carroll S. Ultrabithorax regulates genes at several levels of the wing-patterning hierarchy to shape the development of the Drosophila haltere. Genes & Development. 1998; 12(10):1474–82.

    Article  CAS  Google Scholar 

  61. Cripps RM, Black BL, Zhao B, Lien CL, Schulz RA, Olson EN. The myogenic regulatory gene mef2 is a direct target for transcriptional activation by twist during Drosophila myogenesis. Genes & Development. 1998; 12(3):422–34.

    Article  CAS  Google Scholar 

  62. Held Jr LI. Imaginal Discs: the Genetic and Cellular Logic of Pattern Formation, vol. 39. New-York: Cambridge University Press; 2005.

    Google Scholar 

  63. Abouheif E, Wray GA. Evolution of the gene network underlying wing polyphenism in ants. Science. 2002; 297(5579):249–52.

    Article  CAS  PubMed  Google Scholar 

  64. Tomoyasu Y, Arakane Y, Kramer KJ, Denell RE. Repeated co-options of exoskeleton formation during wing-to-elytron evolution in beetles. Current Biology. 2009; 19(24):2057–065.

    Article  CAS  PubMed  Google Scholar 

  65. Shbailat S, Khila A, Abouheif E. Correlations between spatiotemporal changes in gene expression and apoptosis underlie wing polyphenism in the ant Pheidole morrisi. Evolution and Development. 2010; 12(6):579–90.

    Article  CAS  Google Scholar 

  66. Carroll SB, Gates J, Keys DN, Paddock SW, Panganiban G, Selegue JE, et al. Pattern formation and eyespot determination in butterfly wings. Science (New York, NY). 1994; 265(5168):109.

    Article  CAS  Google Scholar 

  67. Shbailat SJ, Abouheif E. The wing-patterning network in the wingless castes of myrmicine and formicine ant species is a mix of evolutionarily labile and non-labile genes. Journal of Experimental Zoology Part B: Molecular and Developmental Evolution. 2013; 320(2):1–10.

    Article  CAS  Google Scholar 

  68. Peterson T, Müller GB. What is evolutionary novelty? process versus character based definitions. Journal of Experimental Zoology. Part B: Molecular and Developmental Evolution. 2013; 320(6):345–350.

    Google Scholar 

  69. Conway Morris S. Life’s Solution: Inevitable Humans in a Lonely Universe. New-York: Cambridge Univ Press; 2003.

    Book  Google Scholar 

  70. Gould SJ. Wonderful Life: the Burgess Shale and the Nature of History: Random House; 2000.

  71. Colosimo PF, Hosemann KE, Balabhadra S, Villarreal G, Dickson M, Grimwood J, et al. Widespread parallel evolution in sticklebacks by repeated fixation of ectodysplasin alleles. Science. 2005; 307.5717:1928–33.

    Article  CAS  Google Scholar 

  72. Jeong S, Rokas A, Carroll SB. Regulation of body pigmentation by the abdominal-b hox protein and its gain and loss in drosophila evolution. Cell. 2006; 125(7):1387–99.

    Article  CAS  PubMed  Google Scholar 

  73. Protas ME, Hersey C, Kochanek D, Zhou Y, Wilkens H, Jeffery WR, et al. Genetic analysis of cavefish reveals molecular convergence in the evolution of albinism. Nature genetics. 2005; 38(1):107–11.

    Article  PubMed  CAS  Google Scholar 

  74. Sucena E, Delon I, Jones I, Payre F, Stern DL. Regulatory evolution of shavenbaby/ovo underlies multiple cases of morphological parallelism. Nature. 2003; 424(6951):935–8.

    Article  CAS  PubMed  Google Scholar 

  75. Tenaillon O, Rodríguez-Verdugo A, Gaut RL, McDonald P, Bennett AF, Long AD, et al. The molecular diversity of adaptive convergence. Science. 2012; 335(6067):457–61.

    Article  CAS  PubMed  Google Scholar 

  76. Green DA, Extavour CG. Convergent evolution of a reproductive trait through distinct developmental mechanisms in drosophila. Developmental biology. 2012; 372(1):120–30.

    Article  CAS  PubMed  Google Scholar 

  77. Manceau M, Domingues VS, Linnen CR, Rosenblum EB, Hoekstra HE. Convergence in pigmentation at multiple levels: mutations, genes and function. Philosophical Transactions of the Royal Society B: Biological Sciences. 2010; 365(1552):2439–450.

    Article  CAS  Google Scholar 

  78. Travisano M, Mongold JA, Bennett AF, Lenski RE. Experimental tests of the roles of adaptation, chance, and history in evolution. Science. 1995; 267(5194):87–90.

    Article  CAS  PubMed  Google Scholar 

  79. MacLean RC, Bell G. Divergent evolution during an experimental adaptive radiation. Proceedings of the Royal Society of London Series B-Biological Sciences. 2003; 270(1524):1645–50.

    Article  Google Scholar 

  80. Blount ZD, Borland CZ, Lenski RE. Historical contingency and the evolution of a key innovation in an experimental population of escherichia coli. Proceedings of the National Academy of Sciences. 2008; 105(23):7899.

    Article  CAS  Google Scholar 

  81. McGhee GR. Convergent Evolution: Limited Forms Most Beautiful: MIT Press; 2011.

  82. Wittkopp PJ, Stewart EE, Arnold LL, Neidert AH, Haerum BK, Thompson EM, et al. Intraspecific polymorphism to interspecific divergence: genetics of pigmentation in drosophila. Science. 2009; 326(5952):540–4.

    Article  CAS  PubMed  Google Scholar 

  83. Albertson RC, Streelman JT, Kocher TD, Yelick PC. Integration and evolution of the cichlid mandible: the molecular basis of alternate feeding strategies. Proceedings of the National Academy of Sciences of the United States of America. 2005; 102(45):16287.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  84. Kolbe J, Leal M, Schoener T, Spiller D, Losos J. Founder effects persist despite adaptive differentiation: A field experiment with lizards. Science. 2012; 335(6072):1086–9.

    Article  CAS  PubMed  Google Scholar 

  85. Steiner CC, Römpler H, Boettger LM, Schöneberg T, Hoekstra HE. The genetic basis of phenotypic convergence in beach mice: similar pigment patterns but different genes. Molecular Biology and Evolution. 2009; 26(1):35–45.

    Article  CAS  PubMed  Google Scholar 

  86. Brakefield PM, Gates J, Keys D, Kesbeke F, Wijngaarden PJ, Monteiro A, et al. Development, plasticity and evolution of butterfly eyespot patterns. Nature. 1996; 384(6606):236–42.

    Article  CAS  PubMed  Google Scholar 

  87. Jablonski D. Scale and hierarchy in macroevolution. Palaeontology. 2007; 50(1):87–109.

    Article  Google Scholar 

  88. Charlesworth B, Lande R, Slatkin M. A neo-darwinian commentary on macroevolution. Evolution. 1982;:474–98.

  89. Abouheif E. Parallelism as the pattern and process of mesoevolution. Evolution and Development. 2008; 10(1):3–5.

    Article  PubMed  Google Scholar 

  90. Arthur W. Micro-, macro-, and megaevolution In: Hall BK, Olsen WM, editors. Keywords and Concepts in Evolutional Developmental Biology. Harvard University Press: 2007. p. 496.

  91. Leroi AM. The scale independence of evolution. Evolution & development. 2000; 2(2):67–77.

    Article  CAS  Google Scholar 

  92. Meyer JR, Dobias DT, Weitz JS, Barrick JE, Quick RT, Lenski RE. Repeatability and contingency in the evolution of a key innovation in phage lambda. Science. 2012; 335(6067):428–32.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  93. Erwin DH. Macroevolution is more than repeated rounds of microevolution. Evolution & Development. 2000; 2(2):78–84.

    Article  CAS  Google Scholar 

  94. Gibson G, Hogness DS. Effect of polymorphism in the drosophila regulatory gene ultrabithorax on homeotic stability. Science. 1996; 271(5246):200–3.

    Article  CAS  PubMed  Google Scholar 

  95. Zimmer C, Emlen DJ. Evolution: Making Sense of Life: Roberts; 2013.

  96. Hall B. Strickberger’s Evolution: Jones & Bartlett Learning; 2008.

  97. Ridley M. Evolution, 3rd Edition: Blackwell Publishing; 2003.

  98. Stern DL. Evolution, Development, & the Predictable Genome. Greenwood Village: Roberts and Company Publishers; 2011.

    Google Scholar 

  99. Erwin DH, Davidson EH. The evolution of hierarchical gene regulatory networks. Nature Reviews Genetics. 2009; 10(2):141–8.

    Article  CAS  PubMed  Google Scholar 

  100. Davidson EH, Erwin DH. Evolutionary innovation and stability in animal gene networks. Journal of Experimental Zoology Part B: Molecular and Developmental Evolution. 2010; 314(3):182–6.

    Google Scholar 

  101. Parmesan C. Ecological and evolutionary responses to recent climate change. Annual Review of Ecology, Evolution, and Systematics. 2006;:637–69.

  102. Seager R, Ting M, Held I, Kushnir Y, Lu J, Vecchi G, et al. Model projections of an imminent transition to a more arid climate in southwestern North America. Science. 2007; 316(5828):1181–4.

    Article  CAS  PubMed  Google Scholar 

  103. Brown JL, Li SH, Bhagabati N. Long-term trend toward earlier breeding in an american bird: A response to global warming?. Proceedings of the National Academy of Sciences. 1999; 96(10):5565–569.

    Article  CAS  Google Scholar 

  104. Inouye DW, Barr B, Armitage KB, Inouye BD. Climate change is affecting altitudinal migrants and hibernating species. Proceedings of the National Academy of Sciences. 2000; 97(4):1630–33.

    Article  CAS  Google Scholar 

  105. Bhatkar A, Whitcomb WH. Artificial diet for rearing various species of ants. Fla. Entomol. 1970; 53:229–32.

    Article  Google Scholar 

  106. Moreau C, Bell C, Vila R, Archibald S, Pierce N. Phylogeny of the ants: Diversification in the age of angiosperms. Science. 2006; 312.5770:101–104.

    Article  CAS  Google Scholar 

  107. Zhu Y, Bergland AO, González J, Petrov DA. Empirical validation of pooled whole genome population re-sequencing in Drosophila melanogaster. PloS one. 2012; 7(7):41901.

    Article  CAS  Google Scholar 

  108. Skalski GT, Couch CR, Garber AF, Weir BS, Sullivan CV. Evaluation of dna pooling for the estimation of microsatellite allele frequencies: a case study using striped bass (Morone saxatilis). Genetics. 2006; 173(2):863–75.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  109. Vos P, Hogers R, Bleeker M, Reijans M, Vandelee T, Hornes M, et al. Aflp - a new technique for dna-fingerprinting. Nucleic Acids Research. 1995; 23(21):4407–414.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  110. Bonin A, Ehrich D, Manel S. Statistical analysis of amplified fragment length polymorphism data: a toolbox for molecular ecologists and evolutionists. Molecular Ecology. 2007; 16(18):3737–758.

    Article  CAS  PubMed  Google Scholar 

  111. Zhivotovsky LA. Estimating population structure in diploids with multilocus dominant dna markers. Molecular Ecology. 1999; 8(6):907–13.

    Article  CAS  PubMed  Google Scholar 

  112. Vekemans X, Beauwens T, Lemaire M, Roldán-Ruiz I. Data from amplified fragment length polymorphism (aflp) markers show indication of size homoplasy and of a relationship between degree of homoplasy and fragment size. Molecular ecology. 2002; 11(1):139–51.

    Article  CAS  PubMed  Google Scholar 

  113. Mantel N. Detection of disease clustering and a generalized regression approach. Cancer Research. 1967; 27:209–20.

    CAS  PubMed  Google Scholar 

  114. Beerli P. Comparison of bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics. 2006; 22.3:341–345.

    Article  CAS  Google Scholar 

  115. Beerli P, Palczewski M. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics. 2010; 185(1):313–463.

    Article  PubMed Central  PubMed  Google Scholar 

  116. Haag-Liautard C, Coffey N, Houle D, Lynch M, Charlesworth B, Keightley PD. Direct estimation of the mitochondrial dna mutation rate in Drosophila melanogaster. PLoS biology. 2008; 6(8):204.

    Article  CAS  Google Scholar 

  117. Obbard DJ, Welch JJ, Kim KW, Jiggins FM. Quantifying adaptive evolution in the drosophila immune system. PLoS genetics. 2009; 5(10):1000698.

    Article  CAS  Google Scholar 

  118. Cornuet JM, Ravigne V, Estoup A. Inference on population history and model checking using dna sequence and microsatellite data with the software diyabc (v1.0). Bmc Bioinformatics. 2010; 11(401).

  119. Beaumont MA, Zhang WY, Balding DJ. Approximate bayesian computation in population genetics. Genetics. 2002; 162(4):2025–035.

    PubMed Central  PubMed  Google Scholar 

  120. Felsenstein J. Phylip-phylogeny inference package (version 3.2). Distributed by the author: Department of Genome Sciences, University of Washington, Seattle; 2004.

  121. Correa C, Hendry AP. Invasive salmonids and lake order interact in the decline of puye grande Galaxias platei in western patagonia lakes. Ecological Applications. 2012; 22(3):828–42.

    Article  PubMed  Google Scholar 

  122. Wheeler DE, Nijhout HF. Imaginal wing discs in larvae of the soldier caste of Pheidole bicarinata vinelandica forel (hymenoptera: Formicidae). International Journal of Insect Morphology and Embryology. 1981; 10(2):131–9.

    Article  Google Scholar 

  123. Patel N. Imaging neuronal subsets and other cell-types in whole-mount drosophila embryos and larvae using antibody probes. In: Methods in Cell Biology. San Diego: Academic Press Inc.: 1994. p. 445–87.

    Google Scholar 

  124. Warton DI, Wright IJ, Falster DS, Westoby M. Bivariate line-fitting methods for allometry. Biological Reviews. 2006; 81(2):259–91.

    Article  PubMed  Google Scholar 

  125. Tautz D, Pfeifle C. A non-radioactive insitu hybridization method for the localization of specific rnas in drosophila embryos reveals translational control of the segmentation gene hunchback. Chromosoma. 1989; 98(2):81–5.

    Article  CAS  PubMed  Google Scholar 

  126. Metscher BD. Microct for comparative morphology: simple staining methods allow high-contrast 3d imaging of diverse non-mineralized animal tissues. BMC physiology. 2009; 9(1):11.

    Article  PubMed Central  PubMed  Google Scholar 

  127. Mani MS, Giddings LE. Ecology of Highlands: Dr. W. Junk.; 1980.

  128. Hodkinson ID. Terrestrial insects along elevation gradients: species and community responses to altitude. Biological Reviews. 2005; 80(3):489–513.

    Article  PubMed  Google Scholar 

  129. Harrison RG. Dispersal polymorphisms in insects. Annual Review of Ecology and Systematics. 1980;:95–118.

  130. Peeters C, Molet M. Colonial reproduction and life histories In: Lach L, Parr C, Abbott K, editors. Ant Ecology. New-York: Oxford University Press: 2010. p. 159–76.

    Google Scholar 

  131. Murcia C. Edge effects in fragmented forests: implications for conservation. Trends in Ecology & Evolution. 1995; 10(2):58–62.

    Article  CAS  Google Scholar 

Download references


We thank J. Losos, G. Wray, R. Barrett, T. Sanger, Y. Idaghdour, M.B. Dijkstra, A.P. Hendry, L. Nilson, J-P. Lessard, A. Gonzalez, J. Marcus, K Chaudhary, and Abouheif Lab members for comments. We also thank E. Furlong for the anti-mef2 antibody and R. Keller for technical advice. This work was supported by grants from Fonds de la Recherche du Québec - Nature et Technologies, Canada Research Chairs program, and National Sciences and Engineering Research Council (NSERC) to E.A., NSERC postgraduate scholarship to M-J.F.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Ehab Abouheif.

Additional information

Competing interests

The authors declare that they have no competing interests.

Author’s contributions

MJF designed the study, performed field sampling and laboratory work, collected the data, performed the analyses and wrote the manuscript. RAJ provided samples and field assistance. SC provided field assistance and discovered the species. SH, BDM and GBM performed the CT-scan analyses. GBM contributed substantially to revisions. SG performed genetic analyses. EA designed the study, performed field sampling and wrote the manuscript. All authors read and approved the final manuscript.

Additional file

Additional file 1

Supplementary figures and tables. This files includes supplementary tables and figures to the main manuscript. (PDF 4085kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Favé, MJ., Johnson, R.A., Cover, S. et al. Past climate change on Sky Islands drives novelty in a core developmental gene network and its phenotype. BMC Evol Biol 15, 183 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: