Mountains as Islands: Species Delimitation and Evolutionary History of the Ant-Loving Beetle Genus Panabachia (Coleoptera, Staphylinidae) from the Northern Andes

The ant-loving beetle genus Panabachia Park 1942 is a poorly studied beetle lineage from the new world tropics. We recently collected Panabachia from several previously unrecorded locations in the páramo biome of the high Ecuadorian Andes, with males exhibiting great morphological variation in the distribution of the foveae and depressions in the pronotum, as well as aspects of the male genitalia. Here, we employ phylogenetic and species delimitation methods with mitochondrial (COI) and nuclear protein-coding (wingless) gene sequences to examine the concordance of morphological characters and geography with hypothesized species boundaries. Three methods of species delimitation (bPTP, GMYC and Stacey) were used to estimate the number of species, and divergence times between putative species using molecular clock calibration. Phylogenetic analysis revealed two parallel radiations, and species delimitation analyses suggest there are between 17 and 22 putative species. Based on clade support and concordance across species delimitation methods we hypothesize 17 distinct clusters, with allopatric speciation consistent with most geographic patterns. Additionally, a widespread species appears to be present in northern páramo sites, and some sister species sympatry may indicate other diversification processes have operated on certain lineages of Panabachia. Divergence time estimates suggest that Panabachia originated in the Miocene, but most species analyzed diverged during the Pliocene and Pleistocene (5.3–0.11 Mya), contemporaneous with the evolution of páramo plant species.


Introduction
The Andes mountain chain along the South American spine has a dynamic geological and climatological history. A wide range of geological processes, such as plate subduction, volcanism, crustal shortening and terrain accretion, has shaped the topography and the distribution of the species in the tropical Andes [1][2][3][4][5]. The orogenic formation of the Andes started during the Mesozoic and peaked with a massive uplift over the past 30 Ma [6]. The increase in elevation affected the climatic patterns of the region, leading to Quaternary formation of glaciers on the mountain summits [7,8].
These geological and paleoclimatical events have influenced tropical alpine ecosystems, as well as the distributions and genetic diversity of multiple evolutionary lineages that inhabit the tropical Andes, leading to high numbers of endemic species [3,9,10]. During interglacial fluctuations many suitable habitats moved from the mountain slopes into the inter-Andean valleys [5,[11][12][13], sometimes

Field Collection
Samples for this study were obtained from leaf litter samples from 7 sites across the highlands of Ecuador (Figure 1, Table 1). Three leaf litter samples were extracted per site, from a variety of litter types (Polylepis forest, moss, shrubs and grass). The selection of sites was based on conservation status, since most of the collecting took place within the network of protected areas. Collecting permits for this study were previously obtained (MAE-DNG-ARGG-CM-2014-004). The sifted material was transported to the lab and processed using Berlese funnels into 100% ethanol. Collected beetles were separated into morphospecies, based on characters examined ( Figures S1 and S2). Wing size was also recorded for each specimen to understand their flight ability. Voucher specimens of this study will be deposited in the Museo de Zoología de la Pontificia Universidad Católica del Ecuador (QCAZ) after the study is concluded.

Field Collection
Samples for this study were obtained from leaf litter samples from 7 sites across the highlands of Ecuador (Figure 1, Table 1). Three leaf litter samples were extracted per site, from a variety of litter types (Polylepis forest, moss, shrubs and grass). The selection of sites was based on conservation status, since most of the collecting took place within the network of protected areas. Collecting permits for this study were previously obtained (MAE-DNG-ARGG-CM-2014-004). The sifted material was transported to the lab and processed using Berlese funnels into 100% ethanol. Collected beetles were separated into morphospecies, based on characters examined ( Figures S1 and S2). Wing size was also recorded for each specimen to understand their flight ability. Voucher specimens of this study will be deposited in the Museo de Zoología de la Pontificia Universidad Católica del Ecuador (QCAZ) after the study is concluded.   The entire body of each beetle was used to extract genomic DNA using the GeneJet Genomic DNA Purification Kit (Thermo Fisher Scientific, Vilnius, Lithuania). Polymerase chain reaction was used to amplify two molecular markers: COI and wingless. The mitochondrial gene COI was amplified using the primers C1-J-2183 (5 -CAACATTTATTTTGATTTTTTGG-3 ) and TL2-N-3014 (5 -TCCAATGCACTAATCTGCCATATTA-3 , [39]) following the amplification profile described by Caterino and Tishechkin (2014) [40]. For the nuclear gene wingless, we used the primers wg550f (5 -ATGCGTCAGGARTGYAARTGYCAYGGYATGTC-3 ) and wgAbRZ (5 -CACTTNACY TCRCARCACCARTG-3 , [41]) following the amplification profile described by Parker and Grimaldi (2014) [42]. PCR reactions of 25 µL generally contained 2-3 µL genomic DNA, 17.5 µL water, 2.5 µL 10× buffer, 0.5 µL dNTPs, 0.75 µL MgCl 2 , 0.1 µL AmpliTaq ® DNA polymerase (Thermo Fisher Scientific, Waltham, MA, USA) and 1 µL of each primer (10 µm). Amplification cycles were performed in a Mastercycler ® nexus (Eppendorf). PCR products were purified using ExoSAP-IT (USB/Affymetric, Santa Clara, CA, USA), and sequencing was done commercially by Macrogen USA, Inc. (Rockville, MD, USA). Sequences were manually verified and trimmed using Geneious R8 (Biomatters Ltd., Auckland, New Zealand), and aligned using MAFFT v.7 (http://mafft.cbrc.jp/alignment/server/).

Phylogenetic Analyses
Models of molecular evolution were assessed using JModeltest 2.0 [43] for each molecular marker, where the GTR+I+G model appeared to be in the 100% confidence interval for COI and wingless data. The method of Templeton, Crandall and Sing (TCS) v1.21 [44] was used to construct haplotype networks for each data set. DnaSP V6.12 [45] was used to phase the wingless data set. To reconstruct phylogenetic relationships among haplotypes RAxML version 8.2.8 [46] was launched from Mesquite's Zephyr package [47], with 1000 bootstrap replicates. For Bayesian inference, we used Mr. Bayes 3.2 [48] through 20 million generations under default settings.
Divergence time estimates were generated in BEAST 2.0 [49] using a partitioned two-gene data set. Two points of calibration were used, all designating minimum node age within the Pselaphinae: the first was an undescribed Bythinini from Burmese amber (99 Mya) [42]; the second was Parker and Grimaldi's (2014) estimate that higher Pselaphinae arose 150 Mya. We employed an uncorrelated relaxed clock, using a log-normal distribution, and ran the analysis for 20 million generations. Output trees were generated using TreeAnnotator 2.0.02 (http://beast.bio.ed.ac.uk), with maximum clade credibility (MCC) after a 10% burn-in.

Species Delimitation Analyses
The number of species within the genus Panabachia is unknown. Morphology-based species classification is a useful tool to determine species, yet the number of species can be masked by lack of species-diagnostic characters in one sex. Frequently in Pselaphinae, females do not exhibit diagnostic characters. Such is the case in Panabachia [36,50], and a high proportion of females were found among our samples collected. Therefore, we used three sequence-based species delimitation methods to hypothesize the number of reproductively isolated clades in Panabachia from páramo. The methods employed included: the Bayesian implementation of the Poisson tree processes (PTP) model, using the bPTP server (http://species.h--its.org/ptp/) [51]; a single threshold Generalized Mixed Yule Coalescent (GMYC) using the GMYC server (http://species.h--its.org/gmyc/; [52]; and a multi-species coalescent method, Species Tree And Classification Estimation Yarely (STACEY) v.1.2.4 [53], implemented in BEAST 2.0 [49]. Species delimitation analyses using these three models were performed using single locus and multilocus data sets.
For analyses in bPTP, trees generated in Mr. Bayes were used as input. Analyses were run through 100,000 MCMC generations, with a thinning of 100 and 10% burn-in. For GMYC, ultrametric trees were produced in BEAST 2.0 [49], using an uncorrelated relaxed clock, a constant coalescent speciation process prior, through 10,000,000 generations and 10% burn-in. Effective Sample Size (ESS) was evaluated in Tracer v1.5 [54], considering runs with ESS values above 200. Output trees were generated in TreeAnnotator 2.0.02 (http://beast.bio.ed.ac.uk), using maximum clade credibility (MCC) after a 10% burn-in and median heights for node heights. Resulting trees were used as input in the GMYC server, using a single threshold. Lastly, for the implementation of STACEY in BEAST 2.0 [49], files were generated in BEAUTI v.2.4.0 [49]. For the minimal number of clusters, two scenarios were analyzed from all taxa's specimens divided as different species to clusters defined by site. The epsilon value was set to 1 × 10 −4 , following guidelines in the software documentation; nucleotide substitution models were estimated a priori using PAUP 4.0 [55]; a fossilized birth-death model was selected for speciation; and the uncorrelated lognormal model was used to describe the relaxed molecular clock. Input files were run for 500 million iterations, sampling every 10,000th generation. Two replicates were run for each set, ESS values were evaluated in Tracer v.1.6, and independent runs were combined using LogCombiner v.2.4.0 [49] after a 10% burn-in; output trees were summarized in TreeAnnotator [49].

Sequence Data and Polymorphisms
Sampling from litter resulted in the collection of 68 adult Panabachia from seven localities ( Figure 1). For phylogenetic analyses, we used 10 individuals per site where possible (Table 1), although Releche only yielded seven individuals. The COI gene was amplified from 67 samples (GenBank accessions MN536369 to MN536434, Table 2), and the alignment of this gene had a total of 765 base pairs. Of the 765 base pairs, 240 were variable, 66 were parsimony informative (Table 3), and 30 distinct haplotypes were identified using TCS ( Figure 2). The wingless gene was amplified from 62 individuals (GenBank accessions MK674898 to MK674959, Table 2). The alignment for this gene had 445 base pairs with 108 segregating sites from which 105 were parsimony informative (Table 3). A total of 65 alleles were identified in the wingless data set ( Figure 3). Phased data showed that most individuals are homozygotic, with only 16 heterozygotic individuals. Amongst the two data sets, only a few haplotypes/alleles were shared among sites. Individuals from La Virgen shared COI haplotypes with individuals from Pichincha (H9), Atillo (H11) and Releche (H10) (see Figure 2). For the wingless data set, three alleles are shared among sites. One allele is shared by the northern populations of La Virgen, Cayambe and Pichincha (H10), and La Virgen shares alleles with el Angel (H11) and Cayambe (H15; Figure 3). The overall nucleotide diversity for each data set was low (wingless: π = 0.047, COI: π = 0.113; Table 3).

Species Delimitation Analyses
Results from single locus and multilocus analyses using three models of species delimitation are summarized in Figure 4. Multilocus analyses identified 17-22 putative species, with a high level of congruence for most groups. The single locus analyses, however, showed a high variation among outputs. Wingless showed a wide range of results depending on the method used: STACEY suggested 62 clades (out of 68 individuals); bPTP showed 51 clades but GMYC suggested only three species-level clusters ( Figure S3). For the mitochondrial gene, results from bPTP and STACEY showed a higher level of correspondence with results from the multilocus analyses, with 20 clades identified in bPTP and 22 in STACEY. This was not seen in the GMYC analysis, where only four clades were observed ( Figure S1).  For further interpretation of phylogenetic relationships and geographical distribution of genetic clusters, a conservative number of putative species (17) will be considered the most reasonable hypothesis (Figure 4). This hypothesis was based in part on the statistical support for each lineage in the phylogenetic inferences and the distribution of haplotypes in the TCS network analysis. We also assume that some splits observed in bPTP and STACEY multilocus analyses probably represent intraspecific variation; for example, for clade 11 from Mojanda (bPTP) and clade 3 from Cayambe (STACEY). This suggestion is based on the low number of mutations among these haplotypes in TCS analyses. This number of presumed species also corresponds well with preliminary morphological analyses: characters in the pronotum of the male, such distribution of the foveae and shape of the pronotum, show correspondence with the proposed groups ( Figures S1 and S2). As of yet, not all putative species are represented by male specimens, so complete correspondence cannot be assessed. We were also able to trace the presence and absence of wings to assess potential dispersal ability for each proposed clade (Figure 4). We found that most clades are represented by macropterous specimens (clade 6-12 and 16-17). However, numerous wing polymorphic clades are found (clade 1-3, 5 and 13-15), where males are macropterous and females are micropterous. One clade was represented only by micropterous individuals (clade 4).

Phylogenetic Analyses and Divergence Time Estimates
Phylogenetic inferences generated through Maximum Likelihood and Bayesian methods show two separate clades in the combined data set tree as well for the COI gene tree (Figures 4 and 5). Most putative species identified through species delimitation analyses appear to be well supported by bootstrap and posterior probability values (Figure 4). Exceptions are found for clades 1, 2 and 17, which do not have strong branch support (Figure 4). Gene trees did not recover all the genetic clusters observed in the species delimitation analyses (Figures 5 and 6). For example, the COI gene tree shows strong bootstrap and posterior probability values for most putative species ( Figure 5). But clades 2 and 3 are exceptions. In the case of the wingless gene tree, bootstrap and posterior probability values only support seven clades out of the 17 proposed (Figure 6), and multiple incongruences were found in this gene, contributing to the lower phylogenetic resolution, where clades 2, 3, 4, 5, 7, 11, and 14 appear at multiple points. In contrast, the COI tree only presented one incongruence, for individuals of clade 2. The COI tree topology largely resembles the combined data set result, where clade 2 is resolved as a single clade (Figure 4). In reference to the distribution of putative species, clear geographical splits among sites were not detected in the phylogeny, since multiple species were identified for the majority of sites (e.g., four clades in Mojanda, Figure 4). Furthermore, the two main clades both include representatives across most of the area sampled ( Figure 4). However, at a finer scale, some sister clades do show interesting allopatric disjunctions. This was the case of clade (Pichincha, W) and 7 (Atillo, E), separated by distance and side of the mountain range. The separation between clades 9 (Mojanda, W) and 10 (El Angel, W), separated by <70 km, spans the dry Mira River valley.
Divergence time estimates suggest that Panabachia in Ecuadorean páramo originated in the Miocene (9.2 Mya, Figure 7), and that most of the proposed species diverged during the Pliocene and Pleistocene (5.3-0.11 Mya). Two parallel radiations are observed through phylogenetic analyses; these radiations started in the Miocene (8-7.85 Mya) and continued throughout the Pleistocene. The first radiation event (7.86-0.47 Mya) gave rise to clade 7-17, and its sister gave rise to clade 1-6 (4.65-0.24 Mya). This second radiation is composed mostly of northern groups (clades 1-4, and 6), with the exception of clade 5, represented by specimens collected in Atillo (SE). Clade 2 was also exceptional, in that it contained individuals from several sites (La Virgen, Pichincha, El Angel and Releche), the last of which is quite remote from the others. In reference to the distribution of putative species, clear geographical splits among sites were not detected in the phylogeny, since multiple species were identified for the majority of sites (e.g., four clades in Mojanda, Figure 4). Furthermore, the two main clades both include representatives across most of the area sampled ( Figure 4). However, at a finer scale, some sister clades do show interesting allopatric disjunctions. This was the case of clade (Pichincha, W) and 7 (Atillo, E), separated by distance and side of the mountain range. The separation between clades 9 (Mojanda, W) and 10 (El Angel, W), separated by <70 km, spans the dry Mira River valley.
Divergence time estimates suggest that Panabachia in Ecuadorean páramo originated in the Miocene (9.2 Mya, Figure 7), and that most of the proposed species diverged during the Pliocene and Pleistocene (5.3-0.11 Mya). Two parallel radiations are observed through phylogenetic analyses; these radiations started in the Miocene (8-7.85 Mya) and continued throughout the Pleistocene. The first radiation event (7.86-0.47 Mya) gave rise to clade 7-17, and its sister gave rise to clade 1-6 (4.65-0.24 Mya). This second radiation is composed mostly of northern groups (clades 1-4, and 6), with the exception of clade 5, represented by specimens collected in Atillo (SE). Clade 2 was also exceptional, in that it contained individuals from several sites (La Virgen, Pichincha, El Angel and Releche), the last of which is quite remote from the others.

Discussion
High elevation species are particularly interesting given the climatic diversity, high levels of isolation and complex geological history of mountain systems [4,9,56,57]. Yet, alpine beetle faunas from the Andes have only been superficially explored. Previous work in the Ecuadorian páramo has shown distinct patterns of genetic distribution in ground beetles, from higher population structure in flightless ground beetles [28] to high levels of genetic connectivity between populations of a macropterous species [27]. Still, the genetic diversity of other alpine beetle lineages from the Andes has not been assessed, leaving the question open as to whether other alpine insect lineages are following similar patterns as the ground beetles.
Species delimitation analyses facilitate the identification of distinct evolutionary lineages within a sample of individuals [51,52,58], and the use of multilocus genetic data has proven to be a powerful tool for delimiting species [58]. Yet, methods to delimit species vary greatly in parameters and

Discussion
High elevation species are particularly interesting given the climatic diversity, high levels of isolation and complex geological history of mountain systems [4,9,56,57]. Yet, alpine beetle faunas from the Andes have only been superficially explored. Previous work in the Ecuadorian páramo has shown distinct patterns of genetic distribution in ground beetles, from higher population structure in flightless ground beetles [28] to high levels of genetic connectivity between populations of a macropterous species [27]. Still, the genetic diversity of other alpine beetle lineages from the Andes has not been assessed, leaving the question open as to whether other alpine insect lineages are following similar patterns as the ground beetles.
Species delimitation analyses facilitate the identification of distinct evolutionary lineages within a sample of individuals [51,52,58], and the use of multilocus genetic data has proven to be a powerful tool for delimiting species [58]. Yet, methods to delimit species vary greatly in parameters and outcomes, and the search for congruence across results from species delimitation models can provide more reliable hypotheses for species boundaries. Results from the species delimitation analyses show that Panabachia from páramo comprises a diversity of species; 17-22 putative species were identified with bPTP, STACEY and GMYC using a multilocus data set. The three models of species delimitation showed similar outcomes for most species clades (Figure 4). Clustering disagreement was only reported in four clades, where bPTP and STACEY tended to subdivide lineages more finely, showing 1-5 more clades than GMYC. Evidence from phylogenetic inferences suggests that some delimited species using bPTP and STACEY might actually represent intrapopulation variation. This seems to be the case for clade 11 from Mojanda, which is divided by bPTP into two species lineages. In clade 5, from Atillo, bPTP also recognizes two presumed species. This appears be the result of the occurrence of highly divergent genotypes in the wingless data set (3 haplotypes in the wingless data set vs 1 haplotype in the COI data set). For clades 2 and 3, a total of one, two, or four species clades are recognized, depending on algorithm (bPTP, GMYC, or STACEY, respectively). Based on clade support and distribution of haplotypes in TCS, these individuals appear to represent two distinct species clusters.
More broadly, páramo Panabachia appear to represent two parallel radiations, at least by COI and the combined data (Figures 4 and 5), though relationships within each are not clearly resolved. Divergence time estimates show that these radiations started during the Miocene (5.59-7.81 Mya), but 14 out of 17 putative species of Panabachia from páramo originated during the Pleistocene (0.11-4.6 Mya, Figure 7). These estimates are contemporary with the environment they live in, since the Andes reached its current elevation during the Pleistocene [1]. The increase in elevation created suitable conditions for the development of high elevation species [59], which are thought to have evolved from closely related lineages from the lowland tropical areas, as well as from lineages from temperate regions [9,60,61]. Studies of plant lineages from páramo also show accelerated rates of diversification during this period of time [9].
The main factors that affect cladogenetic events are associated with geographical isolation or ecological shifts [62]. Although most of the genetic clusters of Panabachia are well supported by bootstrap and posterior probability values, and represent distinct geographical areas, few phylogenetic relationships between clades are adequately supported to connect to geological events or features. Some correspondence with geography was found in the COI gene tree ( Figure 5). For example, the split between species 8 and 7 spans opposite sides of the mountain range (Atillo and Pichincha). However, most of the speciation events with well supported branches show divergence between presumed species on the same sides of the mountain ranges but separated by distance (clades 9-10 and 14-15). In some instances, for example between clade 9 and 10, the distance might be reinforced by the presence of a putative barrier (Mira river and valley).
While allopatric speciation might explain some of the patterns in Panabachia from páramo, it does not explain the high number of sympatric clades found in El Angel (4), Mojanda (4) and Atillo (3). Climatic oscillation during the Pleistocene glaciations and topographic characteristics of each cordillera appear to influence the level of connectivity and fragmentation of populations across the northern Andes [5]. Understanding the timing of local geological events might give us insight into the factors that influence diversity in these sites; for example, volcanic activity dates to 10,000 years BP in Mt. Mojanda and Mt. Fuya in Mojanda, and Mt. Chiles in El Angel [63,64]. For Atillo, evidence for the presence of small glaciers during the Pleistocene and Holocene was found on Mt. Ayapungo (4730 m) and Mt. Coyay (4630 m) [65], mountains that are adjacent to Atillo. The increased number of putative species of these particular sites could be the result of multiple re-colonization events from adjacent areas during recent environmental fluctuations, as seen in other mountain systems [66,67].
Apart from the effect of environmental conditions and ecological interactions, the dispersal ability of each beetle lineage plays an important role in species diversification [28,[68][69][70]. From previous studies done in ground beetles from páramo, we understand that the loss or reduction of wings can promote diversification events for some beetle lineages [28]. This appears not to be a factor for most putative species of Panabachia from páramo, where more than half of the proposed species are macropterous (56%, Figure 4). Still, wing polymorphic species of Panabachia represent a substantial proportion of the assessed clades (41%, Figure 4), where males are macropterous and females are micropterous or brachypterous. Yet, distinct wing morphologies are not restricted to a geographical area, and both wing polymorphic and macropterous groups can be found in the same sites. Further sampling and a better understanding of the distributional range and natural history of each species is needed to assess whether dispersal ability is a significant driver of diversification of this beetle lineage.
The high diversity found in Panabachia might be associated with the high diversity of leaf litter types sampled. These included decomposing grass leaves and roots, Polylepis and Compositae leaf litter, and moss over rocks and rotten wood. This microhabitat diversity is related to the high diversity of plants and plants forms (from cushion plants and shrubs, to herbaceous rosettes) páramo presents [60,71]. However, microhabitat-focused sampling has not been systematic enough to measure correspondence of putative species of Panabachia with particular types of leaf litter. More focused sampling and more information about the natural history of the group will be needed to determine if any such associations are related to diversification rates.
In comparison with patterns found among ground beetles from páramo, Panabachia from páramo diverged in more recent times (mostly in the Pleistocene). Most ground beetle species from Ecuadorian páramo evolved during the Miocene, prior to the evolution of this ecosystem (~6-20 Mya) [27,28]. Phylogeographic breaks in Panabachia are not as clear as in the Dercylus lineage (Carabidae, Harpalinae), where the presence of geographical barriers (e.g., rivers, dry valleys and mountain range) had a great effect on the pattern of speciation of this group [28]. Yet, most of the proposed species within the Panabachia lineage do represent restricted geographical areas.
When the diversity of Panabachia is compared to the widely distributed Pelmatellus columbianus (11.19 Mya, Carabidae, Harpalinae), clade 2 of Panabachia show similar patterns in the distribution of the genetic diversity, since members of this clade are present across multiple sites. Preliminary analyses of the population structure of species 2 (not shown), suggest Releche represents a distinct genetic cluster. We found incongruence among gene trees for an individual from Releche (SIMT284), which in the COI gene tree is grouped with clade 2, while in the wingless tree is found as sister to clade14. Hence, a more comprehensive analysis will require additional samples from each population to determine if this clade represents one or two proposed species. Nevertheless, widely distributed species of ground beetles have been reported at the same sites as this widespread clade of Panabachia. Such is the case in Pelmatellus columbianus (Cayambe, La Virgen, Pichincha, Releche) and Dyscolus alpinus (Cayambe, Pichincha, La Virgen [28], which suggests similar factors are affecting the distribution of these northern beetle lineages. In particular, the effect of Quaternary glaciation might have enabled gene flow between sites now isolated by elevation [14]. The study of multiple beetle lineages from páramo is slowly providing a better understanding or the evolution of high elevation beetle faunas. An increasing number of studies have concentrated on ground beetle species from the Ecuadorian Andes [17,72,73]. Other beetle lineages from páramo have been less studied, but results to date show that high elevation faunas tend to have an elevated number of endemic species [14,17,[72][73][74]. Numerous putative species of Panabachia (those documented herein as well as others) have yet to be described, further underscoring the importance of conserving high elevation ecosystems. Although most of the sampled sites are already protected areas [75], many high elevation areas across the Ecuadorian Andes are not part of this network of national parks.

Conclusions
Overall, Panabachia represents a promising model for the study of diversification of beetles from high elevation areas in the Andes, considering its high interpopulational diversity and recent divergence (Pliocene and Pleistocene, 5.3-0.11 Mya). Through the use of three methods of species delimitation analyses and two molecular markers we were able to identify 17 putative species from seven sites in the Ecuadorian Andes. Assessed species appear to be restricted to small geographical ranges, with exception of one species clade present in multiple sites in the northern Ecuadorian Andes. The distribution of genetic diversity of Panabachia is complex, and a generalized pattern for alpine beetles of the Ecuadorian Andes has yet to emerge. Multiple factors appear to be shaping the genetic diversity of this beetle lineage, such as mountain isolation, habitat discontinuity and dispersal capability. The pattern of divergence observed across the tree topologies in this study certainly does not capture the entire genetic diversity of the Panabachia lineage, since sampling was only focused on isolated páramo patches. Further sampling (especially cloud and montane forest) and more information about the natural history of the species are needed to develop a more comprehensive picture of the distribution of phylogenetic diversity of Panabachia. The present study should provide a strong preliminary framework for more thorough systematic treatment of this diverse genus.