Congruence between morphology-based species and Barcode Index Numbers (BINs) in Neotropical Eumaeini (Lycaenidae)

Morphology-based species identifications

The basis for identifying the species analyzed in this study is the checklist of Robbins (2004), which includes 1,058 species of Eumaeini in 83 genera. The checklist was updated using subsequent publications (e.g. Balint & Faynel, 2008; Prieto et al., 2008; Duarte & Robbins, 2010; Faynel et al., 2011; Prieto, 2011; Faynel, Busby & Robbins, 2012; Robbins, Heredia & Busby, 2015; Prieto et al., 2016; Prieto & Vargas, 2016; Busby et al., 2017). When necessary, identifications were verified through genitalic examination.

Sampling and sequencing

Collecting permits in Colombia were obtained from ANLA Agencia Nacional de Licencias Ambientales (00594 April 26th 2018). Tissue samples were taken from pinned Eumaeini (Theclinae) in the research collections of Carlos Prieto (RCCP) and Christophe Faynel (RCCF). We selected specimens collected in the past 10 years because older material is more likely to have degraded DNA. Samples came from Costa Rica, French Guiana, Colombia, Ecuador, Peru, and Brazil (Fig. 1). One to three legs were removed from each sampled specimen. The sample included 2,214 specimens of 541 species identified a priori based on the existing classification. The number of specimens per species ranged from 2 to 23.

DNA extraction, amplification, and sequencing of the COI barcode region were carried out by the Canadian Centre for DNA Barcoding (CCDB), Ontario, Canada, using standard high throughput protocols (Ivanova, DeWaard & Hebert, 2006; DeWaard et al., 2008). PCR amplification with a single pair of standard primers targeted a 658 bp region near the 5′ terminus of the mitochondrial cytochrome c oxidase I (COI) gene that included the standard 648 bp barcode region for the animal kingdom (Hebert et al., 2004a). Complete specimen data including images, voucher deposition, accession numbers, GPS coordinates, sequence and trace files are accessible in the Barcode of Life Data System (BOLD dataset: DS-CPCF Faynel-Prieto Neotropical Theclinae; doi: https://dx.doi.org/10.5883/DS-CPCF). Distance-based Neighbor joining (NJ), available on the BOLD website, was used to construct DNA barcode gene trees and to quantify sequence divergence. We analyzed the entire dataset and each genus with the NJ algorithm. In some cases, nearest neighbor genera with few species were combined in a single tree.

Congruence between morphology and BINs

BOLD currently contains close to 9,000,000 barcodes and over 700,000 BINs generated with the Refined Single Linkage (RESL) algorithm. RESL employs a three phased analysis to reach decisions on the number and circumscription of BINs (= MOTUs) in the sequence data set on BOLD (Ratnasingham & Hebert, 2013). It is much faster than other approaches, such as the generalized mixed Yule-coalescent model (Pons et al., 2006; Fujisawa & Barraclough, 2013), a critical requirement for the analysis of large data sets.

Morphological species were partitioned into three categories following the comparative methodology of Hausmann et al. (2013): (I) those in which there was a perfect match between morphological species and BINs; (II) splits: those were morphological species placed in more than one BIN and (III) merges: those where different species shared the same BIN assignment or mixtures where some individuals of a species shared a BIN with another morphological species. We re-examined each sample in the latter two cases by checking both the morphological identification and the alignment and trace files.

Barcode gaps

We analyzed barcode gaps to evaluate the hypothesis that incongruence between morphological species and BINs increase in high Andean lycaenid genera (e.g. Prieto, Nuñez & Hausmann, 2018; Faynel, 2019). The “barcode gap” is a comparison of intraspecific versus interspecific divergence among barcode sequences. A barcode gap exists if the intraspecific divergence (of a particular species) is smaller than its lowest interspecific divergence. For example, a small intraspecific divergence combined with a large interspecific divergence is a large gap. We compared these divergences in the entire dataset and in groups of genera partitioned by the elevational distribution of its species. The criterion for assigning elevational groups was that at least 90% of the species in a genus occur in (1) high mountain habitats (+2,200 m), (2) middle mountain habitats (1,220–2,200 m), (3) lowland habitats (0–1,200 m), (4) middle mountain + high mountain and, (5) or lowland + middle mountain.

As a quick visualization of barcode gaps, we made scatterplots showing maximum intraspecific variation plotted against the minimum distance to the nearest non-conspecific individual. A 1:1 relationship is the point at which the difference between the two is zero (Collins & Cruickshank, 2013). To determine sampling size bias, we also made scatterplots with the number of individuals in each species plotted against their maximum intraspecific variation. These analyses were performed for genera with at least two species and five sequences and for the groups of genera according to their elevational category. To evaluate if the divergence patterns for intraspecific variation and distances to nearest neighbor differ between the sets of species occurring in high mountain and lowland habitats, a Shapiro–Wilk normality test (Shapiro & Wilk, 1965), and a Kruskall–Wallis test were performed. The analyzes were carried out in R software.

Barcode gaps

The percentage of species with a barcode gap in the complete data set was 87.2%. However, the proportion of species from habitats at different elevations in the datasets affected barcode gap frequency. Gaps were observed in 95.7% of the species from lowland ecosystems (0–1,200 m) while only 61.7% of species from exclusively high montane ecosystems (>2,200 m) had clear barcode gaps (Figs. 2, 3). The trend towards a higher percentage of barcode gaps in lowland species was also found when including genera with species distributed in both lowlands and mid-montane habitats (89.5%), exclusively mid montane habitats (82.8%) and genera with exclusively mid- and high-montane species (73.8%) (Table 2). The divergence values to the nearest neighbor were significantly lower in the high Andean species (H = 26.6, p-value = 0.0000002436), while the intra-specific divergence values were significantly lower in the lowland species (H = 6,49, p-value = 0.01084).

BIN sharing

In this study, we obtained 1,839 sequences of 512 species distributed in 90 genera for 557 BIN numbers, representing 78% of the available data on BOLD for the Eumaeini. From the set of specimens that yielded COI fragment lengths of at least 650 bp, 75% of the a priori morphology-based identifications were unambiguously assigned to a single Barcode Index Number (BIN). After a taxonomic a posteriori review, the percentage of perfect matching rose to 85%. Very low levels of interspecific barcode variation can reflect overlooked synonymy if misidentifications are ruled out (e.g. Puillandre et al., 2011), but low genetic divergence, particularly based on just one genetic locus, does not automatically invalidate established taxonomy. In cases of recent phylogenetic divergence, phenotypic differentiation can occur more rapidly than the complete sorting of mtDNA into the new, separated lineages. The decision to consider two species names as synonyms must be made by a taxonomist. That is why in our study we strived for accurate identification before and after delimiting the species using molecular data. When species pairs with low barcode divergences are recovered as monophyletic groups in the cladograms or identification trees, and their morphology is highly divergent, they can be validated as two different species, particularly if they are sympatric. There is no fixed threshold level of divergence that indicates species status because the percentage of divergence that would indicate whether two entities belong to the same species depends on the taxonomic group being studied and its evolutionary history. Nevertheless, most studies have found that COI divergences rarely exceed 2% within named and morphologically validated species, while members of different species typically show higher divergences, and it has been shown repeatedly that this ‘threshold’ can be used in many or most metazoans to determine species status (Ratnasingham & Hebert, 2007, 2013). However, distance and geographic isolation are two aspects that must be taken into account when delimiting biological entities based on established thresholds. Two entities living in sympatry can be considered different species even when there are small genetic divergences of 2% or less. But if these same entities are geographically distant, a 2% divergence may be considered irrelevant to define them as separate species.

The percentage of clearly different morphological species grouped within the same BIN was predominantly high for the montane genus Rhamma, where the BIN BOLD:ABX0547 is shared by five well-differentiated morphospecies, most of them flying in sympatry. We suggest that most of the cases of BIN sharing between morphologically divergent high mountain species represent recently separated lineages that are still undergoing genetic differentiation. As most of these cases were recovered as monophyletic clades in the trees, a lower basic threshold setting in the algorithm would separate these species into different BIN numbers. However, it should be noted that with such modified settings and parameters the number of cases of BIN discordance in the same group of species may increase. In the case of the genus Rhamma the assignment of a single BIN number for clearly different morphological species is a result of the basic settings chosen for delimitation of sequences into BINs.

Incomplete lineage sorting is relatively common in recently and rapidly radiating species groups as these species often have not yet had the necessary time to fix alternative haplotypes or alleles (Galtier & Daubin, 2008). As a result, the relationships of incipient species typically progress from initial polyphyly through paraphyly and reach monophyly once lineage sorting is complete in the two sister species. Thus, in mtDNA analyses, relatively young species may appear polyphyletic or paraphyletic owing to incomplete lineage sorting (Tang et al., 2012). This phenomenon seems to be particularly common in high Andean genera such as Rhamma and has important effects on species identification and delimitation based on genetic sequence analysis. Further studies comparing genetic distances of sympatric and allopatric populations of several pairs of species can help to detect evidence of incomplete lineage sorting, and its prevalence in high mountain species of Theclinae.

BIN splits

High levels of intraspecific barcode variation often reflect cryptic species (e.g. Puillandre et al., 2011, 2012). However, deep barcode splits can also be the result of the recovery of pseudogenes, as a consequence of hybridization, or Wolbachia infection (Werren, Baldo & Clark, 2008; Huemer et al., 2018; Mally, Huemer & Nuss, 2018). High percentages of BIN splits were found in some genera with typical mid- and high mountain representatives such as Podanotum, Johnsonita, Thaeides and Rhamma.

As noted by Prieto, Nuñez & Hausmann (2018), the genus Rhamma includes several species presenting both types of discordance, BIN sharing and BIN splitting (e.g., R. arria and R. bilix). Species with a wider geographical distribution in high Andean ecosystems, seem to show a greater number of incongruences. Morphologically identified specimens were placed in three well-differentiated BINs for R. arria (BOLD:ABX0547, BOLD:ADD3784, BOLD:ADD3785) and four BINs for R. bilix (BOLD:ACF3699, BOLD:ABX0491, BOLD:ADD1839, BOLD:ABX0493). These clades might correspond to either divergent conspecific lineages, or unconfirmed putative species separated, in some cases, by deep genetic divergences. Cases of mitochondrial introgression can hinder the delimitation of some Eumaeini species in the genus Calycopis (Cong et al., 2017), and we suppose that such processes occur more frequently in high Andean genera. Introgressive hybridization may have been common throughout the evolutionary history of these genera which are, therefore, of particular interest to taxonomists and evolutionary biologists because partial and unequal gene exchange can have important effects on the dynamics of speciation and phylogenetic patterns (Grant, 1998; Grant, Grant & Petren, 2005; Funk & Omland, 2003), and affect species identification and delimitation based on DNA sequences.

Barcode gap analysis

A useful display of distance data for species delimitation is a scatterplot showing for each species the maximum intraspecific variation against the minimum distance to the nearest non-conspecific species (‘nearest neighbor’), with a 1:1 slope representing the point at which the difference between the two is zero (Collins & Cruickshank, 2013). This type of representation shows the barcode gap for each species in the dataset and can be an accurate display of the percentage of species in the study group that have a barcode gap (Figs. 2, 3). Since the identification of species and the delimitation of taxonomic entities using barcodes, depend on the existence of a clear DNA barcode gap, a quick visualization of the existence of these gaps in each species can indicate the usefulness of the DNA barcoding approach in a given genus.

Major topographic and climatic variations are important factors that determine the nature of South American biodiversity. The geography of South America, together with its climate and biodiversity, have evolved over a very long period, initiated about 100 MY ago on the ancient Gondwanan continent, although the Andes rose much later (about 25 MY ago). Many species of butterflies are found near, or within, the complex valley systems of the Andes, a result of the combination of at least two important factors: altitudinal gradients and geographic barriers created by the intricate system of valleys and ridges (e.g. Holzinger & Holzinger, 1994; Willmott, 2003; Ebel et al., 2015). The tectonic rise of the Andes has created new environments and modified others, and the uplift of the cordilleras has separated butterfly communities favoring the evolution of allopatric vicariants. Dramatic changes in global climate during glaciations, accompanied by major adjustments in vegetation, created new biomes which again may have stimulated the evolution of new species and subspecies, specially at high altitudes (Purser, 2015; Pyrcz et al., 2017). These changes are quite rapid on a geological scale, and certain subspecies, notably among the high altitude pronophilines (Satyrinae), seem to have evolved since the last glacial maximum 20,000 years ago (Adams, 1985; Pyrcz, Wojtusiak & Garlacz, 2009; Casner & Pyrcz, 2010). Several phylogenetic studies of butterflies indicate that the most recent diversification events tend to occur at high elevations and that the highest altitude species and subspecies are the youngest (e.g., Jiggins et al., 2006; Casner & Pyrcz, 2010; Pyrcz et al., 2017).

Due to incomplete lineage sorting in very young species, it is not easy to accurately define taxonomic boundaries, and additional difficulties may be caused by hybridization. Although the RESL algorithm as the basis for the BIN system (Ratnasingham & Hebert, 2013) provides a powerful tool to propose primary, tentative species hypotheses for large datasets (Ratnasingham & Hebert, 2007), such an mtDNA-based approach cannot prove the absence of gene flow and still depends on arbitrary, a priori settings and assumptions. The efficiency of these methods largely depends on the accumulation of mtDNA mutations since species separation, and thus can only delimit lineages with sufficiently long isolation (Rannala, 2015). Moreover, nothing is known about how the kind of speciation process (vicariance of an existing species versus a small founder colonization) might affect the ability of barcodes to identify species correctly. We assume that incomplete lineage sorting and occasional hybridization are usual phenomena in the very young species of South American Eumaeini, and that these are the two most likely causes of the low percentages of DNA barcode gaps found in high Andean species in comparison with the older species from the lowlands.

However, other very species-rich groups of Andean butterflies with recent speciation processes such as the subtribe Pronophilina, have shown very high percentages of barcode gaps and perfect congruence between morphology and DNA barcodes (e.g. Marín et al., 2017; Pyrcz et al., 2018). This shows that, in certain groups, other biological factors allow young high Andean species to present more complete DNA lineage sorting in short periods of time. In the case of Pronophilina, the low vagility of its species could be a determining factor that limits gene flow between separate populations and promotes lineage sorting.