The role of dispersal for shaping phylogeographical structure of flightless beetles from the Andes

Study area

The study area encompassed the Ecuadorian Andes, where all the sampling took place above 3,100 m in the páramo ecosystem. The selection of sites was based on previous records of occurrence of the target taxa, and potential geographical barriers. Some of these potential barriers include: dry valleys, rivers, faults, and the east/west split of the Andes mountain range (Krabbe, 2008; Guayasamin et al., 2010; Quintana et al., 2017). In particular, the Pastaza, Chanchan, and Mira rivers, and their corresponding valleys, appear to limit the distributions of some bird, amphibian, and plant species (Krabbe, 2008; Guayasamin et al., 2010; Quintana et al., 2017). The Pallatanga fault, a prominent NE–SW strike–slip fault crossing the western cordillera (Baize et al., 2015), may also affect lineage dispersal. The most prominent geographic divide considered in this study was the split of the Ecuadorian Andes into east and west cordilleras, associated with limits of distribution in several vertebrate and plant species (Chaves et al., 2007; Guayasamin et al., 2010; Sklenář, Dušková & Balslev, 2011).

Sample collection

Ground beetle specimens were collected at eight sites along the Ecuadorian Andes, between 3,100 and 4,000 m in national parks during two field seasons in the summers of 2015 and 2016 (Fig. 1; Table 1). Permissions necessary for the collection of samples were previously obtained through the Ministry of Environment of Ecuador (MAE), permit number MAE–DNG–ARGG–CM–2014–004). Dyscolus alpinus was found in six out of the eight sites sampled (Fig. 1A), while Dercylus specimens were collected in five sites (Fig. 1B). Adult beetles were collected into 100% EtOH using three methods: manual collection (underneath rocks and on vegetation), pitfall trapping, using two 100 m transects (one trap every 10 m), where traps ran for 2–3 days; and through the collection of leaf litter samples from the páramo floor. Three samples of leaf litter were collected per site, and were processed in the lab using Berlese funnels. Ground beetles were identified using Moret’s (2005) taxonomic keys. An additional sample was provided and identified by Pierre Moret as Dercylus cordicollis (SIMT265) from Cotacachi, this sample was included for the analysis of Dercylus cordicollis populations. Outgroup taxa were chosen based on availability of fresh material and taxa considered closely related carabid lineage within the Harpalinae. Outgroup taxon sampling included 12 taxa: two were collected in the Ecuadorian páramo; six taxa were collected in the southeastern United States; and four additional samples were retrieved from GenBank (KR604910, Table S1). Species in the genera Platynus, Agonum, and Incagonum were used as outgroups in the analyses of Dyscolus. Amara, Oodes, Dicaelus, Sarticus, Notonomus, and Loxodactylus were used as outgroups in the analyses of Dercylus (Table S1). Voucher specimens will be deposited in the Museo de Zoología de la Pontificia Universidad Católica del Ecuador (QCAZ) after the study is completed.

DNA extraction, sequencing, and alignment

Genomic DNA was extracted from 110 specimens (with an average of 10 individuals per site, per species) using GeneJet’s Genomic DNA Purification Kit (Thermo Fisher Scientific, Vilnius, Lithuania) following the manufacturer’s instructions. A fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene was amplified using the primers TL2–N–3014 (5′–TCCAATGCACTAATCTGCCATATTA–3′ and C1–J–2183 (5′–CAACATTTATTTTGATTTTTTGG–3′; Simon et al., 1994) and the amplification profile described in Caterino & Tishechkin (2015). One nuclear protein-coding gene, CAD, was amplified using the primers CD439F (5′–TTCAGTGTACARTTYCAYCCHGARCAYAC–3′) and CD688R (5′–TGTATACCTAGAGGATCDACRTTYTCCATRTTRCA–3′). Heminested PCRs were done for difficult samples, using CD439F (5′–TTCAGTGTACARTTYCAYCCHGARCAYAC–3′) and CD1098R (5′–TTNGGNAGYTGNCCNCCCAT–3) following Wild & Maddison (2008). All PCRs were performed in 25 μl reactions, containing 17.5 μl water, 2.5 μl 10x buffer, 0.5 μl dNTPs, 0.75 μl MgCl2, 0.1 μl DreamTaq DNA Polymerase (Thermo Fisher Scientific, Vilnius, Lithuania) and one μl of each primer (five nm). PCR products were cleaned using ExoSAP–IT (USB/Affymetric, Santa Clara, CA, USA). Sequencing was done commercially by Macrogen USA, Inc. (Rockville, MD, USA).

DNA sequences were edited using Geneious R8 (Biomatters Ltd., Auckland, New Zealand) and aligned in MAFFT v.7 (Katoh & Standley, 2013) together with sequences from outgroups. Models of sequence evolution were tested using JModeltest 2.0 (Darriba et al., 2012) for each data set, where GTR + I + G model was in the 100% confidence interval for both COI and CAD data.

Genetic diversity and population structure

The levels of genetic diversity were calculated using Arlequin 3.5.2.2 (Excoffier & Lischer, 2010) and DnaSP (Rozas et al., 2003) using the default settings. Groups analyzed were: (1) Dyscolus alpinus; (2) Dercylus orbiculatus; and (3) Dercylus species together (Dercylus orbiculatus¸ Dercylus cordicollis, Dercylus praepilatus and Dercylus sp. from Atillo and Cotacachi) in order to contrast inter– and intrapopulation variability. Parameters estimated included: estimated haplotype diversity (h), nucleotide diversity (π), number of polymorphic sites (S), number of migrants per generation (Nm), and Tajima’s D neutrality test statistic (D). ɸ_ST values were calculated to assess the population structure, while an analysis of the molecular variance (AMOVA) was performed to test the influence of geographical barriers, including the rivers, dry valleys discussed above, and the division between populations of eastern and western cordilleras. A Mantel test of correlation between genetic and geographic distances was performed using 1,000 randomizations (also in Arlequin). Lastly, Structure 2.3.4. (Pritchard, Stephens & Donnelly, 2000) was used to estimate the number of populations (K), using no admixture for the mtDNA and admixture for the nuclear gene. A total of 10 iterations were run per potential K, and the presence of one to eight populations were tested for both beetle lineages.

Phylogenetic analyses and the molecular clock

Phylogenetic analyses were performed to determine relationships among individuals for each lineage and among lineages. To explore marker-specific effects, mitochondrial and nuclear data sets were initially analyzed separately, but combined data were used to generate the most robust hypotheses for each beetle lineage using MrBayes 3.2 (Huelsenbeck & Ronquist, 2001) and RAxML version 8.2.8 (Stamatakis, 2014). The analyses performed in MrBayes incorporated a GTR + G + I model, with two independent runs through 1,000,000 generations, with trees sampled every 100th generation. For RAxML (version 8.2.8; Stamatakis, 2014), analyses were launched from Mesquite’s Zephyr package (Maddison & Maddison, 2015) with 1,000 bootstrap replicates using default settings. In addition to these phylogenetic inferences, the number of haplotypes was estimated using DnaSP (Rozas et al., 2003) and haplotype networks were inferred with The method of Templeton, Crandall and Sing (TCS) (Clement, Posada & Crandall, 2000), using default settings.

To understand the timing and rate of diversification for each ground beetle lineage, divergence times were estimated using an uncorrelated relaxed clock in BEAST 2.0 (Bouckaert et al., 2014) using the COI and CAD data matrixes. For Dyscolus alpinus, one point of fossil calibration was used to designate minimum node ages, using a log-normal distribution and enforced monophyly. This fossil represents the same tribe as Dyscolus (Platynini): a fossil of Limodromus from the late Eocene (37.2–33.9 Mya), of Baltic origin (Schmidt, 2015). The divergence time for the Dercylus lineage was also estimated using two points of fossil calibration: Amara (Harpalinae: Zabrini) from the Middle Eocene (44.1 Mya), of Baltic origin (Larsson, 1978), and Pterostichus walcotti (Harpalinae: Pterostichini) from the Eocene (37.2–33.9 Mya) of Teller County, Colorado (Scudder, 1900). In addition to fossil calibrations, substitution rates presented by Andújar, Serrano & Gámez-Zurita (2012) for COI were used in both analyses for comparison. Log files from BEAST were examined in Tracer 1.5 (http://tree.bio.ed.ac.uk/software/tracer/) to ensure effective samples sizes. Trees were visualized in TreeAnnotator 2.0.02 (http://beast.bio.ed.ac.uk), using maximum clade credibility after a 10% burn–in.

Maximum entropy models for species in the genus Dercylus

The ecological niche of each species of Dercylus was modeled to determine if species distribution overlapped, using present bioclimatic variables at 2.5’ spatial resolution (WorldClim v.1.4; Fick & Hijmans, 2017) with the use of MAXENT 3.3.3 (Phillips, Anderson & Schapire, 2006). Dercylus from Cotacachi and Atillo were excluded from these analyses given the few records available. For other species of Dercylus, sites mentioned in Moret (2005) were georeferenced using GeoLocate (Rios & Bart, 2010) and combined with the distribution data generated in this study. A total of 41 occurrence records were gathered for analyses. The model of distribution was based on five bioclimatic variables: BIO₁ = Annual Mean Temperature, BIO₂ = Mean Diurnal Range (Mean of monthly (max temp–min temp)), BIO₅ = Max Temperature of Warmest Month, BIO₆ = Min Temperature of Coldest Month, BIO₁₂ = Annual Precipitation, BIO₁₃ = Precipitation of Wettest Month and BIO₁₄ = Precipitation of Driest Month. The analyses ran for 10 replicates of 500 iterations and using a 10% training presence, using the default settings for the convergence threshold (10⁻⁵). For each species, the output from the average run was reclassified using ArcMap 10.4.1 (Esri, Redlands, CA, USA) with a 10% threshold. Model performance was evaluated using the area under the curve (AUC) calculated by MAXENT, where values between 0.7 and 0.9 indicated a good discrimination (Swets, 1988).

Genetic diversity and population structure in Dyscolus alpinus

The portion of the COI gene analyzed for Dyscolus alpinus was 767 bp long, with 41 parsimony informative sites making up six haplotypes (Fig. 2A; GenBank accessions MK440253–MK440258). Haplotype 1 is the most widespread across populations, found in four of six sites (Fig. 2A). Unique haplotypes were recorded from Atillo (H3), Mojanda (H5), and Cayambe (H6; Fig. 2A). Haplotype diversity (h) is shown in Table 2 for each population, where Mojanda exhibits the highest haplotypic diversity (h = 0.71). The overall nucleotide diversity of COI was low (π = 0.00–0.01). Structure analyses tended to overestimate the number of populations; results from multiple simulations show that the number of populations using this data set is one (K = 1). Results from the neutrality test for the overall values (D = –0.81, p = 0.42; F = –2.16, p = 0.40), as well as on most populations are not significant (Table 2). Cayambe and Atillo were exceptions, but these did not exhibit a high proportion of rare alleles, often associated with population expansion (Fig. 2A). Additionally, a Mantel test showed no correlation (p = 0.13, R² = 0.08) between genetic (ɸ_ST) and geographical distances.

The CAD data set for Dyscolus alpinus spanned 744 bp, with 64 parsimony informative sites. These represented five haplotypes (GenBank accessions MK440260–MK440264), with no heterozygote individuals observed. Haplotype 1 is found in four of the six populations. Unique haplotypes were recorded from Atillo (H4) and Cayambe (H5). Haplotype networks for both genes show that most haplotypes are present at multiple sites (Fig. 2B). The overall haplotypic (h) and nucleotide diversity was low for CAD data set (Table 2). Among all the sites analyzed, Mojanda showed the highest haplotypic diversity (h = 0.91). Structure analyses of the nuclear data set showed three populations (K = 3). Similar to the results from the mitochondrial data, most neutrality tests using CAD (overall value D = –0.86; p = 0.36; F = –5.97, p = 0.40) were not significant, with Cayambe and Pichincha as exceptions. Again, however, these apparent exceptions exhibit few haplotypes, inconsistent with a population expansion scenario. Results from the Mantel test show there was no correlation between genetic and geographic distances (p = 0.81, R² = 0.01).

The AMOVA of nuclear and mitochondrial data for Dyscolus alpinus shows that most of the variation occurs within populations (Table 3). These results seem to contrast with the overall ɸ_ST values (COI ɸ_ST = 0.54; CAD ɸ_ST = 0.26), which suggest moderate levels of population subdivision among populations. When ɸ_ST values were compared on a population–by–population basis, some genetic differentiation was found among populations (Table 4; Fig. 3). COI data exhibited higher overall genetic differentiation between sites (ɸ_ST = 0.32–0.90; Table 4). Some exceptions are found between sites that are closer in distance (less than 36 km), such as Cayambe–Mojanda (separated by 30 km): ɸ_ST = 0.10, and La Virgen-Pichincha (separated by 36 km): ɸ_ST = –0.10. The CAD data set shows lower genetic differentiation among all populations of Dyscolus alpinus (Table 4). However, higher ɸ_ST values are observed when Cashca Totoras is compared with the other sites (ɸ_ST = 0.5–0.65; Fig. 3), except when compared to Cayambe (ɸ_ST = –0.08). The overall number of migrants between Dyscolus populations (Nm = 0.31) indicates that there is population subdivision across sites. Higher levels of genetic connectivity were observed between Mojanda and Cayambe (Nm = 1.79), results that are also supported by the low ɸ_ST values between these two sites (Table 4).

Genetic diversity and population structure in species of the genus Dercylus

In the case of Dercylus species, the portion of the COI gene analyzed was 776 bp, with 76 parsimony informative sites constituting 16 haplotypes (GenBank accessions MK440220–MK440235). No haplotypes were shared among sites (Fig. 4A), and overall nucleotide diversity within these species is low (π = 0.00–0.01; Table 2). Eight haplotypes were recorded for populations of Dercylus orbiculatus, whereas two haplotypes were recovered from Dercylus cordicollis. The individual sampled from Cotacachi represents a distinct lineage of Dercylus closely related to Dercylus cordicollis, based on the TCS analyses and phylogenetic inference, though only one haplotype was found. For specimens of Dercylus praepilatus from Salinas, one haplotype was recorded, while four haplotypes were recorded from specimens collected from the locality of Atillo. Based on the morphological similarity with Dercylus orbiculatus, specimens from Atillo were thought to be another population of Dercylus orbiculatus; however, based on haplotype designation and other phylogenetic analyses, this population appears to be a distinct lineage/species of Dercylus. Haplotype diversity indices are shown in Table 2 for each species and population, where populations of Dercylus orbiculatus from Culebrillas show the highest haplotypic diversity (h = 0.86) among Dercylus species and populations. Structure analyses of Dercylus orbiculatus individuals all resolved as one population, even though the two sites do not share haplotypes. These analyses were also performed for the entire clade, and one population was found (K = 1). Neutrality test values were not significant (overall value D = 0.26, p = 0.77; F = 0.89, p = 0.44), with the exception of the population of Dercylus orbiculatus from Culebrillas, for which significantly positive D and Fs were found (Table 2). Lastly, the Mantel test performed with ɸ_ST and geographic distances show no significant correlation (R² = 0.05, p = 0.3) between genetic and geographical distances.

The portion of the CAD gene amplified for Dercylus species was 669 bp with 24 parsimony informative sites (GenBank accessions MK440236–MK440251). A total of 13 alleles were identified across the Dercylus lineage and the overall nucleotide diversity was low (π = 0.00–0.02). Four alleles were recovered from Dercylus orbiculatus populations, and three alleles from Dercylus cordicollis, with one shared with Dercylus from Cotacachi (H5; Fig. 4B). Four alleles were identified in the population of Dercylus from Atillo, and three from Dercylus praepilatus. Two alleles were shared between Dercylus praepilatus and Dercylus sp. from Atillo (H8 and H11; Fig. 4B). The highest allele diversity is recorded for the Dercylus sp. from Atillo, followed by Dercylus orbiculatus from El Cajas and Dercylus cordicollis (Table 2). Structure analyses for Dercylus orbiculatus populations alone showed only one population (K = 1), whereas results from the analyses of the whole clade show five distinct populations (K = 5), corresponding to each species. Neutrality tests supported significantly positive Fs for Dercylus cordicollis from Pichincha (Table 2). Results from the Mantel test show no significant correlation (p = 0.8, R² = 0.02) between genetic (ɸ_ST) and geographical distances.

The AMOVA showed that most of the genetic variation was found among groups when testing the effect of geographical barriers (Table 5). The high variation among groups is consistent with the existence of multiple species, since species in this genus do not have widespread distributions, with the exception of Dercylus orbiculatus, occurring across the southern Ecuadorian Andes (Fig. 1; Moret, 2005). The highest percentage of variation in the AMOVA when parsing the data by species was seen across the Chimbo and Chanchan dry valleys and rivers, which overlap with another potential geographical barrier, the Pallatanga depression. These results were also supported by the overall ɸ_ST values, which show high levels of differentiation among species and sites (COI ɸ_ST = 0.89; CAD ɸ_ST = 0.95). A high level of genetic differentiation was also observed when ɸ_ST values were compared between populations and species in COI data set (Table 6). In this data set, the lowest ɸ_ST values were observed between populations of Dercylus orbiculatus from Cajas and Culebrillas (ɸ_ST = 0.50).

In the case of the nuclear coding gene CAD, the overall value of ɸ_ST (ɸ_ST = 0.93) shows high levels of population differentiation. When each population was compared against each other, we found moderate to complete population subdivision between some sites (ɸ_ST =0.77–1.00; Table 6), particularly when comparing against populations of Dercylus orbiculatus from Cajas or Culebrillas. Lower ɸ_ST values were found between the two populations of Dercylus orbiculatus Cajas and Culebrillas (ɸ_ST = 0.10), as well as between Pichincha and Cotacachi (ɸ_ST = 0.00), Pichincha and Salinas (ɸ_ST = 0.00), Salinas and Atillo (ɸ_ST = 0.03), and Salinas and Cotacachi (ɸ_ST = 0.00). The number of migrants per generation was calculated for two populations of Dercylus orbiculatus (Cajas and Culebrillas), which revealed low gene flow (Nm = 0.39).

Phylogenetic analyses and divergence time estimations

Phylogenetic analyses were conducted separately for each gene data set using Bayesian and Maximum Likelihood methods. Analyses of the two fragments for Dyscolus alpinus yields similar resolution (Figs. 5A and 5B), with Dyscolus sp. (SIMT208) as sister to the Dyscolus alpinus clade (Fig. 5A). These results are consistent with the TCS haplotype network, which shows most samples are contained within one network, with the exception of COI haplotype 6 (SIMT208–H6; Fig. 2A). Most haplotypes are present at multiple sites (Fig. 5A), with only two unique haplotypes (Atillo-H3 and Mojanda-H5). Haplotype 1 is found at five of six sampling sites (Mojanda, Cayambe, Pichincha, La Virgen, Cashca Totoras); haplotype 2 and haplotype 4 are each found at two (Atillo, Mojanda and Cayambe, Mojanda, respectively; Fig. 2A). The highest haplotypic diversity was found at Mojanda with four haplotypes (H1, H2, H4, H5), while the lowest was at Cashca Totoras and Pichincha with only the widespread H1 haplotype found at each site. A similar pattern was recovered for CAD. First, CAD also suggested that SIMT208 from Cayambe represented a distinct Dyscolus species (Fig. 5B). Within Dyscolus alpinus, four CAD haplotypes were recorded (Fig. 5B). Most sites yielded haplotype 1 (Mojanda, Cayambe, Pichincha, La Virgen, Cashca Totoras, Fig. 2B), with only one unique haplotype (H4) observed from an individual from Atillo. The Mojanda population showed the highest haplotypic diversity (H1, H2, H3), whereas Pichincha (H1) and Cashca Totoras (H1) each only revealed one widespread haplotype. A relaxed molecular clock was used on the combined matrix for Dyscolus alpinus to calculate divergence time estimates for the clade (Fig. 6). These analyses suggest that the Dyscolus alpinus clade originated during the Miocene 6.09 Mya (1.26–13.64 Mya). The timing of splitting–events for the Dyscolus alpinus clade (14.33 Mya, with a confidence interval of 3.63–27.46 Mya; Fig. 6) is slightly older than most páramo plant species (2–5 Mya; Madriñán, Cortés & Richardson, 2013), but contemporary with the evolution of highland species in the northern Andes mountain chain, which arose during the Miocene (Weir, 2006; Hines, 2008).

For Dercylus, the two genes reveal slightly different patterns (Figs. 7A and 7B). The COI tree suggests that most species and sites are monophyletic (Fig. 7A), with Dercylus orbiculatus the exception, appearing as a paraphyletic grade. The CAD gene tree shows less structure associated with site and species, since multiple alleles are shared among species, for example, between Dercylus praepilatus and Dercylus sp. from Atillo (H8, H11), and between Dercylus cordicollis and Dercylus sp. from Cotacachi (H5; Fig. 4B). Yet, in CAD, Dercylus orbiculatus appears as a distinct clade, supported by bootstrap and posterior probability values (Fig. 7B).

The analysis of the combined matrix (COI and CAD; Fig. 8) exhibits clearer resolution among species and populations of Dercylus, though Dercylus orbiculatus appears again as a paraphyletic grade, as in the COI tree. A relaxed molecular clock based on the combined data set shows that the Dercylus clade originated during the Oligocene 29.81 Mya, with a confidence interval of 12.45–50.38 Mya. Thus, most of the sampled species in this study originated prior to the evolution of páramo, with Dercylus orbiculatus (10.74–13.5 Mya, confidence interval of 1.10–28.58 Mya), Dercylus from Cotacachi (16.39 Mya, confidence interval of 3.60–24.38 Mya) and Dercylus sp. from Atillo (6.04 Mya, originating confidence interval of 1.13–12.90 Mya) during the Miocene (16.4–6.04 Mya). Dercylus cordicollis, on the other hand, originated in the Pliocene (3.68 Mya, confidence interval of 0.00–6.96 Mya), and Dercylus praepilatus (2.38 Mya, confidence interval of 0.01–7.23 Mya) in the Pleistocene. Only these last two diversifications were contemporary with the diversification of most páramo plants (2–5 Mya, Madriñán, Cortés & Richardson, 2013).

Niche modeling for species in the genus Dercylus

The ecological niche model for the three described species of Dercylus showed a high predictive ability under a random model, with AUC values between 0.998–0.999. Populations of this ground beetle lineage are mainly present in the Andean region of Ecuador, with potentially overlapping distributions (Fig. 9). Compared to other species analyzed, the maximum entropy model for Dercylus orbiculatus shows that this species is present at high elevation, but its distribution also extends to lower elevations, especially toward the south. In contrast, Dercylus praepilatus has the smallest potential distributional range. Lastly, niche modeling also revealed that suitable habitats for species of this genus also extend into to the Andean regions of Colombia and Peru. To test if the distribution of the species extends into these areas, further sampling is needed.