Abstract

Up to 9 described species of Junonia butterflies occur in the Americas, but authorities disagree due to species similarities, geographical and seasonal variability, and possible hybridization. In dispute is whether Caribbean Junonia are conspecific with South American species. Cytochrome oxidase I (COI) barcodes, wingless (wg) sequences, and Randomly Amplified Fingerprints (RAF) were studied to reveal Junonia population structure in French Guiana, Guadeloupe, Martinique, and Argentina. Phylogenetic analysis of COI recovered 2 haplotype groups, but most Junonia species can have either haplotype, so COI barcodes are ambiguous. Analysis of nuclear wingless alleles revealed geographic patterns but did not identify Junonia species. Nuclear RAF genotyping distinguished 11 populations of Junonia arranged into 3 clusters. Gene flow occurs within clusters but is limited between clusters. One cluster included all Argentinian samples. Two clusters included samples from French Guiana, Martinique, and Guadeloupe and appear to be divided by larval host plant use (Lamiales versus Scrophulariales). Many Junonia taxa were distributed across populations, possibly reflecting patterns of genetic exchange. We had difficulty distinguishing between the Caribbean forms J. zonalis and J. neildi, but we demonstrate that Caribbean Junonia are genetically distinct from South American J. evarete and J. genoveva, supporting the taxonomic hypothesis that they are heterospecific.

1. Introduction

Buckeye butterflies, genus Junonia (Nymphalidae), are an important model system for experimental research in the Lepidoptera [1, 2]. Junonia species have been widely used to study the evolution and development of butterfly wing colour patterns [29]. Experimental tools to manipulate gene expression developed in Junonia are broadly applicable across the Lepidoptera [1014]. Junonia has also been used in studies of insect endocrinology [1517] and has been an important system for examining the evolution of larval host plant preference and tolerance to host plant toxins [1821].

Junonia butterflies are found throughout the Old and New World tropics. In the Western Hemisphere, forms of Junonia occur from southern Canada to Tierra del Fuego [2224] and have a complicated taxonomic history. In 1775, Cramer [25] identified and described two similar species of Junonia, J. evarete and J. genoveva, from Suriname, a Dutch colony on the north coast of South America. The species were described according to the standards of the time (without designated type specimens) and the descriptions were accompanied by hand-tinted plates that reproduced the colours from the original watercolour drawings of specimens of each form (republished in [26]).

In the 20th century, there was considerable disagreement in the scientific community about whether Cramer’s two species were truly distinct [27, 28] or whether all of the specimens belonged to J. evarete [2933]. This is a result of the geographical [29, 34] and seasonal [35] variability of Junonia and the fact that some Junonia forms closely resemble one another [36, 37]. In addition, different Junonia forms share identical karyotypes () [26, 38] and are capable of hybridization and the production of fertile offspring [3941], further complicating the process of assigning names to Junonia specimens. These features make applying the biological species concept [42], phylogenetic species concept [43], or morphospecies concept [44] very difficult in Junonia. Operationally, we use the isolation species concept that defines species as systems of populations such that genetic exchange between these systems is limited or prevented by one or more reproductive isolating mechanisms [45, 46]. Identifying and understanding the reproductive isolating mechanisms operating in Junonia will be of great importance in clarifying Junonia taxonomy.

Authorities who favoured the two-species hypothesis in Junonia called the larger form J. genoveva and the smaller form J. evarete. In 1985, Turner and Parnell [26], using specimens from Jamaica and Florida, USA, verified the existence of two Junonia species in both regions. However, after consulting Cramer’s hand-tinted plates and comparing them to specimens from Jamaica and Florida, Turner and Parnell [26] switched the names so the larger species was now J. evarete and the smaller species was J. genoveva. Neild [22], using specimens from Venezuela (geographically much closer to the type locality of Suriname), also confirmed the existence of two species. However, Neild [22], unsatisfied with Cramer’s [25] published plates (copies of which differ from one another due to variation among the watercolourists who tinted them and differences in how the plates aged), consulted Cramer’s original watercolours and Junonia specimens from many localities in South America. Using this reference material, Neild [22] reversed Turner and Parnell [26] so that the larger species was again J. genoveva and the smaller species was J. evarete. Neild [22] also designated new types for J. evarete and J. genoveva to facilitate future taxonomic work.

Recently, L. Brévignon and C. Brévignon [23, 47, 48] identified 5 Junonia species (J. evarete, J. genoveva, J. wahlbergi, J. litoralis, and J. divaricata) from French Guiana. This represents the most diverse assemblage of Junonia in the New World. There are two forms of Junonia known from the Caribbean Islands, “zonalis” and “neildi,” which were initially recognized as subspecies of mainland J. evarete and J. genoveva, respectively [49] and later as two distinct species: J. zonalis and J. neildi [23]. In recognizing J. neildi and J. zonalis as distinct species, L. Brévignon and C. Brévignon [23] restricted the use of the species epithets J. evarete and J. genoveva to mainland Central and South American forms. If it were confirmed that the Caribbean forms are actually distinct species with respect to Junonia from the mainland, this would explain some of the widespread difficulty of assigning appropriate taxonomic names to specimens from Florida, Jamaica, and elsewhere in the West Indies. Finally, there are two additional Junonia species, J. coenia from North America and J. vestina from the Andes mountains of South America, for a current total of up to 9 species of New World Junonia.

The first molecular phylogenetic approaches to understanding the relationships among the species discussed here established that the New World fauna appears to be monophyletic and that the various New World forms are indeed in the genus Junonia [50, 51] (some authorities had previously placed these species in the related genus Precis) [52, 53]. Unfortunately, these early studies, which incorporated data from both mitochondrial and nuclear loci, had limited taxon sampling, including data from only 3 New World species [50, 51]. More recent studies of the molecular phylogeny of New World Junonia [23, 54, 55], which have better taxon sampling, have focused entirely on the mitochondrial cytochrome oxidase I (COI) locus, which is widely used as a barcoding locus for animal taxa [56, 57]. Based on mitochondrial haplotype sequences, the relationships among many New World Junonia species are ambiguous and most species are not reciprocally monophyletic [23, 24, 54]. The degree to which recent divergences, retained polymorphisms, and/or hybridization events contribute to these patterns in Junonia is unknown because only mitochondrial markers were considered. What is apparent is that there are two very divergent COI haplotype groups (4% sequence divergence between them) present in New World Junonia: Group A, which predominates in South America and is also present in the Caribbean, and Group B, which predominates in North and Central America but which also occurs in the Caribbean and South America [24, 54]. Sequences belonging to each of these haplotype groups can occur in different individuals of the same species at the same locality [24].

The most successful study to date to distinguish between New World Junonia taxa using molecular markers employed a combination of mitochondrial and nuclear markers to examine populations of Junonia in Buenos Aires, Argentina. Borchers and Marcus [24] used DNA sequences from the nuclear wingless gene and anonymous nuclear loci identified by Randomly Amplified Fingerprinting (RAF) (a technique used to assess genetic diversity within populations [5860] and gene flow between populations [61]) in addition to sequences from the mitochondrial COI gene. They identified 3 distinct populations of Junonia from Buenos Aires: one population with dark-coloured wings referred to as J. evarete flirtea [62] which Borchers and Marcus [24] suggested may correspond to J. wahlbergi and 2 light-coloured populations that correspond to J. genoveva hilaris and either a genetically disparate population of J. genoveva hilaris or an undescribed cryptic Junonia species. However, the relationship of these Argentinian forms with J. evarete and J. genoveva from Suriname and French Guiana is not known, so we refer to them as J. “flirtea” and J. “hilaris.

In the current study, we extend the genetic tools that were employed by Borchers and Marcus [24] to Junonia populations from French Guiana and the French Antilles in order to study the distinctiveness of the named taxa within and between these two localities. This will allow an explicit test of the 7-species taxonomic hypothesis (2 species in the French Antilles plus 5 species in French Guiana) of L. Brévignon and C. Brévignon [23] and also detect possible hybridization events between named forms. By using a common set of markers we will also be able to compare these populations to previously studied Junonia from Argentina [24].

2. Materials and Methods

2.1. Specimens and DNA Preparation

A total of 104 Junonia specimens were collected from the wild as adults, reared from wild-collected larvae, or reared from eggs laid by wild-collected adults and frozen at −20°C (Table 1). DNA was isolated from legs removed from each specimen. Some samples (42 specimens) were prepared by the Canadian Centre for DNA Barcoding at the University of Guelph as previously described [23]. The remaining samples (62 specimens) were processed in our laboratory using the Qiagen DNEasy Blood and Tissue kit (Qiagen, Düsseldorf, Germany) as previously described [24], except that the extractions were performed in a Qiagen QIAcube instrument using the standard instrument protocol for purification of total DNA from animal tissue. Extracted DNA was stored at −20°C.

Only 10 μL aliquots of DNA were available for the 42 Junonia specimens processed at the University of Guelph, which was insufficient for the number of experiments we wished to conduct. To produce additional template, whole genome amplification using Illustra Genomiphi V2 (GE Health Care Life Sciences, Pittsburgh, PA, USA) protocol was performed as follows: 1 μL of DNA template and 9 μL of sample buffer were incubated at 95°C for 3 min, cooled to 4°C, mixed with 9 μL of reaction buffer and 1 μL of enzyme, incubated at 30°C for 90 minutes and then 65°C for 10 minutes, and cooled to 4°C. Deionized distilled water was used as the template for a Genomiphi amplification negative control. Genomiphied samples were stored at −20°C.

2.2. Mitochondrial Cytochrome Oxidase I Protocol

Cytochrome oxidase I (COI) PCR products were generated using a seminested two-step amplification with LCO1490 and Nancy primers followed by a reamplification with LCO1490 and HCO2198 (Table 2) [63, 64]. Quick-Load Taq 2X Mastermix (New England Biolabs, Ipswich, MA, USA) was used in PCR reactions with total volumes of 25 μL. Amplification protocols were run on a BioRad MyCycler or S1000 Thermal Cycler (BioRad, Hercules, California, USA) for these and all other PCR amplifications unless otherwise specified. LCO1490/Nancy PCR reaction conditions were 95°C for 5 minutes; 40 cycles of 95°C for 1 minute, 46°C for 1 minute, 72°C for 1.5 minutes; and a final 5-minute extension at 72°C before being placed on a 4°C hold. LCO1490/HCO2198 PCR reaction conditions were 95°C for 5 minutes; 35 cycles of 94°C for 1 minute, 46°C for 1 minute, 72°C for 1.5 minutes; and a final 5-minute extension at 72°C before being placed on a 4°C hold. PCR reactions were evaluated by gel electrophoresis (1% agarose in TAE buffer, 78 V for 1 hour, visualized with ethidium bromide).

Samples that failed to amplify with LCO1490 and HCO2198 were reamplified with M13-uniminibarF1 (miniCOIF) and M13-uniminibarR1 (miniCOIR) (Table 2) [65]. MiniCOI PCR reaction conditions were 95°C for 2 minutes; 5 cycles of 95°C for 1 minute, 46°C for 1 minute, 72°C for 30 seconds; 35 cycles of 95°C for 1 minute, 53°C for 1 minute, 72°C for 30 seconds; and a final 5-minute extension at 72°C before being placed on a 4°C hold. Further reactions were carried out to obtain overlapping PCR products that could be assembled as contigs to obtain additional sequence data. Additional primers were designed to bind to invariant regions of the Junonia COI gene (miniCOIF2 and miniCOIR2 in one reaction and either miniCOIF3 and HCO2198 or miniCOIF2 and HCO2198 (Table 2) in a second reaction) to selectively amplify required sequences. Reaction conditions for these primers were the same as the miniCOI protocol described above.

2.3. Nuclear Wingless Protocol

Wingless PCR products were generated using lepwg1 and lepwg2 primers (Table 2) [66]. Wingless PCR reaction conditions were 94°C for 5 minutes; 40 cycles of 94°C for 1 minute, 46°C for 1 minute, 72°C for 2 minutes; and a final 10-minute extension at 72°C before being placed on a 4°C hold. While these primers typically work well in Junonia [24], the samples analyzed here failed to produce detectable products, likely due to poor preservation of nuclear DNA. These PCR reactions were used as the template for PCR reamplification with miniwgF and miniwgR (Table 2), which we designed to bracket the most informative interval of the Junonia wingless coding sequence (Table 3). Mini-wingless reaction conditions were 95°C for 5 minutes; 40 cycles of 95°C for 1 minute, 57°C for 1 minute, 72°C for 1 minutes; and a final 5-minute extension at 72°C before being placed on a 4°C hold.

2.4. Sequencing

Correctly sized PCR products were sequenced as previously described [24]. Products were sequenced in both directions, usually with the same primers that generated the products. When the miniwgR primer produced poor quality sequences samples were reamplified with miniwgF and T7-miniwgR and sequenced using T7 primer (Table 2). Sequencing reactions were analyzed on an ABI 3730xl automated sequencer and edited using Sequencher 4.6 software [67]. Sequences were trimmed to the appropriate size (Table 3) and aligned in CLUSTALW [68].

2.5. Randomly Amplified Fingerprinting Protocol

Randomly Amplified Fingerprinting (RAF) was used to gather a large multilocus data set [60]. Amplifications were carried out using single fluorescently labelled primers that act as both forward and reverse primers. A product is produced only if the primers bind in the correct orientation and close enough to one another for amplification. The 3 RAF primers, each covalently bound to a 6-FAM fluorescent molecule (Integrated DNA Technologies, Iowa City, Iowa, USA), used in these amplifications were RP2 (5′-/6-FAM/ATGAAGGGGTT-3′), RP4 (5′-/6-FAM/TGCTGGTTCCC-3′), and RP6 (5′-/6-FAM/TGCTGGTTTCC-3′) [59]. Amplifications were performed in triplicate along with positive and negative (distilled deionized water) controls for a total of 954 RAF amplifications. Reaction volumes of 10 μL were used. Samples were run in a BioRad MyCycler Thermocycler under the following reaction conditions: 95°C for 5 minutes; 30 cycles of 94°C for 30 seconds, 57°C for 1 minute, 56°C for 1 minute, 55°C for 1 minute, 54°C for 1 minute, 53°C for 1 minute; and a final 5-minute extension at 72°C before being placed on a 4°C hold. Reactions were shipped at room temperature to the Biotechnology Core Facility at Western Kentucky University (Bowling Green, Kentucky, USA). 10 μL HiDye formamide and 1 μL RX-500 GeneScan Size Standard (Applied Biosystems, Carlsbad, California, USA) were added to each PCR tube upon receipt. The solution was then vortexed for 1-2 seconds and placed in a microcentrifuge at 13,000 rpm for 30 seconds at room temperature. Samples were placed into individual wells on a sequencing plate and incubated at 95°C for 4 minutes in a thermocycler. Following 3–5 minutes on ice, samples were loaded into an ABI 3130 automated sequencer (Applied Biosystems), which was fitted with a 50 cm capillary filled with Pop-7 sequencing polymer for fragment analysis.

2.6. Mitochondrial Cytochrome Oxidase I Analysis

A subset of the samples in the current study was used in a previous barcoding study performed in another laboratory [23]. To ensure that there was no confusion or contamination of DNA samples during transfer we resequenced COI from 17 samples (of 42 transferred) that had been previously sequenced. In all cases, identical sequences were obtained by our laboratory as previously reported [23]. COI sequence alignments were converted to NEXUS format for phylogenetic analysis using several different reconstruction methods (distance, parsimony, and likelihood) that rely on vastly different assumptions about sequence evolution, each of which recovered essentially the same tree. For the sake of brevity, we will only present the maximum likelihood analysis (HKY model, 10 replicate heuristic searches with random number seeds, tree bisection, and reconnection branch swapping algorithm) [69]. Other previously published Junonia COI sequences were included in this phylogenetic analysis [23, 24, 51, 57, 7074]. We also conducted a maximum likelihood bootstrap analysis of this dataset (500 fast addition replicates, collapsing all notes with frequency less than 50%). The aligned COI FASTA sequences generated by this study along with 22 previously published Argentinian Junonia COI sequences [24] were analyzed using Arlequin 3.5 [75]. We employed an AMOVA analysis with the following settings: 1000 permutations, determining the minimum spanning network (MSN) among haplotypes, computing distance matrix, and pair-wise difference with a gamma value of 0. The minimum spanning tree output from AMOVA was put into HapSTAR-0.7 [76], which displays the haplotype network in graphical form. Since analysis in Arlequin requires all sequences to be of the same length, the analysis was first conducted using all samples that amplified using LCO1490/HCO2198 (Figure 2) and then repeated after trimming all sequences to the length of miniCOIF2/HCO2198 (Figure 3). Additional adjustments to the network were made using Canvas X (ACD Systems, Seattle, Washington, USA) such as scaling the population circles to reflect sample size and adding pie charts to reflect the RAF population assignment or geographical location and species.

2.7. Nuclear Wingless Analysis

For the Junonia species sequenced in this study and Argentinian Junonia wingless sequences from a prior study [24], individuals heterozygous for single nucleotide polymorphisms (SNPs) in the coding sequence were identified using sequencing chromatograms and CLUSTALW alignments. For each polymorphism, the genotype of each individual was entered into PHASE 2.1.1 [77] and analyzed using the default settings. PHASE uses the Markov Chain-Monte Carlo method to group coinherited SNPs in order to determine the most probable wingless alleles present in each individual. The most likely alleles identified in PHASE were assigned to each individual and the data was then entered into GENEPOP 4.0.10 [78]. GENEPOP was used to test for genetic differentiation (Exact G test [79]) by determining if the alleles from each subpopulation were drawn from the same distribution. GENEPOP settings used for testing all populations were a demorisation of 10,000, 10,000 batches and 10,000 iterations per batch. Finally, Structure 2.3.3 [80] was used to analyze the wingless data since, unlike GENEPOP [78], Structure does not require the a priori assignment of individuals to specific subpopulations. Population structure exhibited by wingless alleles was analyzed using Structure 2.3.3 [80] with settings for codominant alleles, a 10,000 step burn-in and 1 million Markov Chain-Monte Carlo Method replicates. Ten replicate structure searches tested each of 15 different population models with 1 to 15 subpopulations among the 88 wingless sequences. The maximum log likelihood () for the 10 replicate searches for each population model was used to calculate the posterior probability () of each population model. Haplotype networks of wingless alleles were constructed in the same manner as COI except that PHASE output identifying the most likely wingless genotypes was formatted for input into Arlequin 3.5 [75].

2.8. Randomly Amplified Fingerprinting Analysis

Fragment analysis sample runs were combined with previously studied Argentinian Junonia [24] and analyzed using GENEMAPPER version 3.7 software (Applied Biosystems). An allelic bin size of 3 base pairs was selected in order to detect polymorphic alleles without introducing excessive noise into the analysis associated with small differences in run time between samples. The resulting GENEMAPPER genotypic classifications were exported to an Excel spreadsheet (Microsoft, Redmond, Washington, USA) for further analysis. Bands that appeared in negative control amplifications (of deionized distilled water with no DNA added) were considered artefacts and removed from further analysis for all samples. Within the 3 replicate RAF fragment runs for each primer from an individual butterfly, allele-calling for the presence or absence of the dominant allele at each RAF locus was based on a majority rule determination (at least 2 of the 3 runs had to show the allele for it to be scored as present). Each locus was coded in binary with 0 indicating the absence of an allele and 1 indicating the presence of the allele. Such binary data was analyzed using Structure 2.3.3 software [80] with the same settings as described previously for the wingless data except that in the case of the RAF data set only dominant alleles could be scored. A total of 50 replicate searches were carried out on each of = 1–15 populations, first including only the samples genotyped in this study (primarily from French Guiana and the Caribbean) and then again including the 22 Argentinian specimens genotyped in a previous study [24].

Allele frequencies for each RAF locus were calculated for each population identified in Structure and formatted for input for the CONTML application of PHYLIP 3.5 [81] as implemented in EMBOSS Explorer [82]. CONTML uses a rigorous maximum likelihood algorithm to estimate phylogenies based on allele frequencies. In this model, all divergence between populations is assumed to be due to genetic drift in the absence of new mutations [83]. CONTML trees were exported in NEXUS format and rendered in EvolView [84] for interpretation. A parallel analysis was conducted in CONTML for the RAF data set and the allele frequency data obtained for COI and wingless.

3. Results

3.1. Mitochondrial Cytochrome Oxidase I Results

New COI DNA sequences generated by this project were deposited in Genbank (67 accessions, numbers KJ469059–KJ469126), with the exception of specimens with only COI minibarcode sequence fragments [65], which were submitted to the DNA Databank of Japan (DDBJ, 5 accessions AB935341–AB935345). Full COI barcode sequences, covering the interval between LCO1490 and HCO2198, were recovered from 65 specimens (17 reported previously [23] and 48 new sequences), 15 of which required assembling 3 sequence contigs to obtain the 658 bp sequence. Partial barcode sequences were obtained from miniCOIF/R sequences assembled into contigs with miniCOIF2/R2 sequences (2 specimens), miniCOIF/R sequences assembled into contigs with miniCOIF3/HCO2198 sequences (2 specimens), and miniCOIF2/HCO2198 sequences alone (16 specimens). Overall, some COI sequence was recovered from 90 of the 104 specimens.

Analysis of the COI sequences produced a maximum likelihood phylogenetic tree (Figure 1). As previously reported [24, 54], there are two distinct mitochondrial haplotype groups in New World Junonia. Haplotype group A is found in South American and Caribbean specimens, while haplotype group B includes many North American, Central American, and Caribbean specimens, as well as some South American specimens. A few forms of Junonia appear to be associated with only one haplotype group (group A: the South American forms J. “flirtea” and J. vestina; group B: the North American forms J. coenia and J. “nigrosuffusa”). All other Junonia species (J. divaricata, J. evarete, J. genoveva, J. “hilaris,” J. litoralis, J. neildi, J. wahlbergi, and J. zonalis) were found to include individuals with haplotypes in both group A and group B. COI coding sequences from the two haplotype groups contain no internal stop codons, no insertions or deletions, and few (and generally conservative) nonsynonymous substitutions and show little evidence of heterozygosity (no double bands in PCR products, very few double peaks in sequencing reads to indicate heterozygous sites within PCR products), suggesting that these are not pseudogenes or nuclear copies of mitochondrial DNAs, but true allelic alternatives.

Junonia litoralis from French Guiana and J. neildi from the Caribbean, both of which feed on black mangrove (Avicennia germinans) as larvae, include individuals with mitochondria from both COI haplotype groups. Junonia genoveva from French Guiana, which has often been considered conspecific with J. neildi from the Caribbean, also includes individuals that carry haplotypes in groups A and B. Junonia evarete from French Guiana have exclusively group A COI haplotypes, but J. zonalis from the Caribbean,which has sometimes been considered conspecific with J. evarete, includes individuals from both haplotype groups. In specimens from Guadeloupe, regardless of species, group A haplotypes are relatively rare (5/15, 33% type A). Specimens from Martinique, which is closer to the South American mainland than Guadeloupe, primarily have group A haplotypes (9/12, 75% type A). Specimens from French Guiana, regardless of species, primarily have group A haplotypes (49/56, 87.5% type A) as do specimens from Argentina (19/22, 86% type A).

Analysis of the haplotype networks produced from COI sequences revealed the same general pattern regardless of whether 86full-length 658 bp barcode sequences (Figure 2) or 102 partial 520 bp barcode sequences (Figure 3) were analyzed. Haplotype groups A and B are clearly delineated with 11 nucleotide changes between the genotypes from groups A (J. zonalis specimen LCB164) and B (J. litoralis specimen LCB320) that are most similar to each other. The two networks differ primarily with respect to how some group A genotypes are connected to LCB320, a specimen of J. litoralis located near the center of haplotype group A. There are also some minor rearrangements among the genotypes in haplotype group B between the two networks. There is a strong geographic signal in the COI haplotype network with specimens from the same place often sharing identical or similar genotypes. Many genotypes are found only in French Guiana, the French Antilles, or Argentina and there is only one group B genotype that can be found in all 3 localities (Figure 3(a)).

The haplotype networks also reveal that almost all J. wahlbergi and almost all J. evarete from French Guiana share a single group A genotype that is rare in other Junonia species. Similarly, the vast majority of J. genoveva from French Guiana possess one of 3 disparate genotypes within COI haplotype group A (Figures 2(a) and 3(a)) that are rare in all other Junonia species. Most of the remaining J. genoveva specimens have sequences that are one nucleotide removed from one of the three haplotype A genotypes or carry a genotype from haplotype group B. For most other species of Junonia there are no abundant COI genotypes that are diagnostic for particular species.

3.2. Nuclear Wingless Results

Short 137 bp wingless sequences were recovered from 66 of the 104 specimens. The newly generated sequences were submitted to DDBJ (accession numbers AB935346–AB935395 and AB936758–AB936773). These sequences were analyzed in combination with 22 wingless sequences from Argentinian Junonia specimens [24]. Following CLUSTALW alignments of the sequences of the mini-wingless PCR products with the Argentinian wingless products, 22 single nucleotide polymorphisms (SNPs) were identified in this highly variable region. These sites were confirmed in the chromatograms by the presence of a double peak, which indicates heterozygosity. Of the 22 SNPs, 21 are binary SNPs and 1 SNP contains 3 alternate nucleotides. Analysis of the SNPs in PHASE 2.1.1 [77] identified a total of 35 wingless alleles. The most probable allelic combinations for each Junonia specimen were identified in PHASE and then assigned to each individual for entry into GENEPOP 4.0.10 [78].

GENEPOP was first used to test for genetic differentiation among all populations. Separating Junonia populations solely by geography () or species () and by both geography and species () suggests significant genetic differentiation and distinct wingless allele distributions associated with these two factors. However, separating individual Junonia specimens using mitochondrial COI haplotype alone to define populations shows no statistically significant distinct distribution of wingless alleles associated with mitochondrial haplotypes (). However, when specimens were categorized by geography and COI haplotype (); species and COI haplotype (); or geography, species, and COI haplotype (), the wingless alleles again appear to be drawn from significantly different distributions among subpopulations, likely due to the extremely strong influence of geography and species on the distribution of wingless alleles.

The wingless data was also analyzed in Structure 2.3.3 [80]. The results of the Structure analysis, testing for 1 to 15 subpopulations (Table 4), showed that the model with the highest posterior probability () based on the wingless allelic data is (all samples belonging to a single population). Haplotype networks of wingless alleles show that one allele (a9) that is common in J. genoveva populations from French Guiana (50%) is much rarer in Junonia from both Martinique and Guadeloupe (15%) (Figure 4). A second Junonia wingless allele (a7) is fairly common in both French Guianan (38%) and French Antillean (54%) populations.

3.3. Randomly Amplified Fingerprinting Results

RAF fragments for the RP2, RP4, and RP6 primers were recovered from all 104 Junonia specimens in this study and 22 Argentinian Junonia from our prior work [24]. RAF produces fragments of several different sizes from amplification with a given primer. 43 RP2 loci, 61 RP4 loci, and 18 RP6 loci were identified for a total of 122 variable RAF loci (Table 5). Structure 2.3.3 software [80] was used to test for 1 to 15 subpopulations among the French Guianan and Caribbean Junonia butterflies. The results of this analysis (Table 6) showed that the model with the highest posterior probability () is (samples belonging to 8 separate populations). Curiously, when this analysis was repeated with the inclusion of Junonia samples from Argentina, the model with the highest posterior probability () was . However, the Structure analysis did not distribute the Argentinian samples among the populations previously established for French Guiana and the Caribbean. Instead, all of the Junonia from Argentina were assigned to 1 population, while the samples from French Guiana and the French Antilles were redistributed among 5 populations. This was very curious because in our prior study the Argentinian Junonia, when analyzed by themselves using Structure, were divided into 3 populations [24]. The ability of Structure to detect population subdivision is reduced when sample sizes are very small, but the software is far more sensitive to insufficient numbers of variable markers [85]. This study employs more variable RAF markers (122 loci) than our earlier study (51 loci), but the additional loci are fixed in Argentinian Junonia [24]. When the Argentinian data set is analyzed in isolation both sets of. RAF loci show the model with the highest posterior probability () is (samples belonging to 3 separate populations) as expected. Populations that were unchanged in composition whether or not Argentinian samples and were included in the analysis are populations 1, 4, and 5 (Figure 5). Structure identifies major discontinuities in population structure most readily. When both major and minor discontinuities exist in the same data set, the minor discontinuities can be missed because Structure employs a heuristic search algorithm that explores the solution space rather than calculating an exact solution [80]. In the absence of major discontinuities, the algorithm more readily identifies minor discontinuities. For populations that fused or divided between analyses with and without Argentinian samples, we reanalyzed each group of affected specimens in Structure separately from all other specimens and we used the population assignments from these separate analyses as the definitive group definitions. Populations defined in this way are indicated by a shared number followed by letters (e.g., 2A and 2B).

Overall, there are 11 populations established by Structure analysis of RAF alleles (Figure 5(a)). RAF Population 1 (yellow) primarily includes specimens of J. litoralis, but 2 J. zonalis, 1 J. genoveva, and 1 J. evarete specimens also show genetic influence from this population.Population 2A (dark blue) includes all specimens of J. wahlbergi and J. divaricata and all but one of the J. evarete specimens. Population 2A also includes individuals from all of the other French Guianan and French Antillean Junonia species. Most of the specimens in Population 2B (light blue) also show genetic influence from Population 2A, but there are 2 Caribbean J. zonalis specimens whose primary genetic influence is population 2B. Population 3A (dark purple) is comprised of J. genoveva and J. zonalis specimens. Populations 3B (light purple) and 3C (pink) both include J. genoveva, J. neildi, and J. zonalis specimens. Population 4 (red) contains only J. genoveva specimens. Population 5 (orange) is primarily composed of J. genoveva specimens, although there is 1 J. zonalis specimen that also shows influence from this population. Populations 6A-C are the same as the 3 Argentinian Junonia populations previously described [24].

Finally, of particular interest is the distribution of specimens showing the genetic influence of more than one population. The populations exist in 3 distinct clusters with some genetic exchange within each cluster, but little apparent genetic exchange between clusters (Figure 5(a)). Populations 1, 3A, 3B, 3C, 4, and 5 belong to one such cluster with Population 3C as the “hub” population, showing genetic exchange with all of the other “satellite” populations, which have varying amounts of (often very limited) genetic exchange with one another. It is interesting to note that the 3C hub population is almost exclusively composed of specimens from the French Antilles. Specimens of J. zonalis and J. neildi from the Caribbean and J. genoveva and J. litoralis from French Guiana comprise the remainder of this cluster (along with 1 specimen of J. evarete). In this cluster, many individuals have RAF genotypes with influence from more than 1 population, suggesting possible past or current genetic exchange among these populations. A second cluster includes Populations 2A and 2B, which show extensive genetic exchange between them, but limited genetic exchange with all of the other populations we have identified. The 2A-2B cluster shows a strong genetic influence from J. evarete, J. wahlbergi, and J. divaricata from French Guiana, although it also includes specimens from the other 2 French Guianan species, specimens of J. neildi and J. zonalis from the Caribbean, and J. coenia from Florida. A third cluster includes Populations 6A, 6B, and 6C and contains only Argentinian specimens. Like the 1-3A-3B-3C-4-5 cluster, one population (6B) has genetic exchange with the other two populations, but there may be little or no direct gene flow between 6A and 6C [24].

CONTML analysis of the RAF allele frequency data for populations identified in Structure produced an arbitrarily rooted tree that shows much of the same clustering that was apparent from other analyses of this data set (Figure 5(b)). Populations 2A and 2B, which contain many Caribbean specimens, are sister groups on the tree while Populations 6A, 6B, and 6C, which consist of only Argentinian specimens, arise as an unresolved polytomy from a single common ancestor. The remaining populations, which make up the 1-3A-3B-3C-4-5 cluster, are also adjacent to one another on the CONTML tree, although they do not form a monophyletic group. CONTML analysis of the RAF allele frequency data combined with COI and wingless allele frequency data produces a tree of identical topology with minuscule changes in branch length (not shown).

4. Discussion

We hoped that using a set of molecular tools that had worked well to distinguish between forms of Junonia from Argentina [24] would allow us to unambiguously distinguish between forms from the French Antilles and French Guiana. Unfortunately, this is not entirely the case. Cytochrome oxidase I (COI) barcodes are clearly unreliable for distinguishing among most New World forms of Junonia because individuals of most named species contain haplotypes from either of the two main haplotype groups (A and B, Figure 1). This is consistent with the findings of prior studies that have found to be identical haplotypes in Junonia [23, 24, 54]. Certain COI haplotypes are found primarily in particular geographic regions, but most are found in more than one species in that region and, thus, are not diagnostic (Figures 2 and 3). There are examples of significant apparently intraspecific mitochondrial haplotype divergence in other species and are typically associated either with the existence of cryptic species (each with different haplotypes) [86, 87] or with hybridization followed by introgression of mitochondria [88, 89]. The presence of cryptic species is usually corroborated by the demonstration of either consistent phenotypic differentiation between cryptospecies, population subdivision based on nuclear markers between cryptospecies, or both. Neither this study nor previous studies of Junonia have been able to provide corroboration supporting the existence of cryptic species associated with COI haplotype diversity [23, 24, 54], so we feel weight of the evidence is more consistent with a history of hybridization and introgression. Except for certain special cases in which hybridization between forms of Junonia is apparently absent from a region [55], COI barcoding should probably be abandoned as a method for distinguishing taxa among New World Junonia. At the same time, the segregation of haplotypes A and B, separated by approximately 4% sequence divergence, in Junonia populations across much of the Western Hemisphere is a phenomenon worth studying in its own right.

There appears to be a north-south gradient across the Caribbean with respect to the frequency of COI haplotypes with the type A haplotype ranging from 0% in North America [54] to 87% in the South American mainland (this study and [24]). There is a long-standing hypothesis that Caribbean Junonia are a ring species [29, 34], which is defined as a group of species or subspecies that exhibit a ring-like distribution such that the forms at the extreme ends of the range overlap [42]. Gene flow occurs through intermediate forms in the middle of the ring, but, in a classic ring species, forms found in the overlapping region of the ring do not interbreed [90, 91]. Western Cuba and the nearby Island of Pines were identified as the possible region of overlap for the putative Junonia ring species since several phenotypically distinct forms of Junonia coexist there [34, 92]. Our results suggest that if there is a region of secondary contact between two ends of a Caribbean Junonia ring species, then hybridization between the terminal forms appears to be occurring and the zone of overlap is not limited to Cuba. More extensive mapping of the distribution of the haplotypes in combination with phylogenetic analysis of additional mitochondrial sequence data may provide insights into the origin of these haplotypes (currently unknown [51]) and what evolutionary forces may be contributing to their continued presence in Junonia populations.

Alleles from the nuclear locus wingless, which were helpful in distinguishing between Junonia from Argentina [24], were less effective in this study. This is probably due, at least in part, to the small 137 bp sequence fragment of wingless that we were able to recover from Caribbean and French Guianan Junonia (versus 402 bp previously recovered for Argentinian Junonia [24]). While the wingless fragment that was recovered is the most variable portion of the New World Junonia wingless coding sequence (Figure 4), potentially informative sequence variation elsewhere in the gene was not available for us to study. It would be highly desirable to obtain additional specimens of Junonia from French Guiana and the Caribbean with better-preserved nuclear DNA so that larger portions of the wingless gene could be analyzed. Another factor that may make wingless sequences less useful in French Guiana and the Caribbean is the large number of forms of Junonia coexisting and possibly interbreeding in this region [23]. If wingless coding sequences or sequences closely linked to wingless are adaptive, for example, contributing to colour pattern phenotypes under selection [8, 93], such sequences may be subject to introgression after hybridization events [94, 95]. As such, the evolutionary history of the introgressed region of the genome may not be representative of the evolutionary history of the organism as a whole [96]. Finally, wingless signalling may contribute to the development of colour pattern phenotypes that are used as field marks for identifying species of Junonia [8], most notably the prominent pale median stripe on the ventral hind wing of some species [97, 98]. This connection between developmental processes and species-specific phenotypes in combination with transgenic techniques [1012, 14] may permit us to identify the specific mutations responsible for phenotypic evolution in Junonia and to characterize the molecular mechanisms responsible for colour pattern diversity.

In the past, Randomly Amplified Fingerprint (RAF) genotyping was extremely effective at distinguishing between forms of Junonia from Argentina [24]. In this study, these populations remained distinguishable from one another (Populations 6A, 6B, and 6C) and were genetically distinct from populations sampled from other regions (Figure 5). In fact, the Argentinian populations are all more genetically similar to one another than to any populations in French Guiana or the French Antilles. Based on wing colour patterns, Argentinian Junonia have been conventionally referred to as J. genoveva hilaris (C. & R. Felder), the light buckeye butterfly, and J. evarete flirtea (Fabricius), the dark buckeye butterfly [62]. However, Borchers and Marcus [24] identified 3 genetically distinct Argentinian populations based on RAF genotypes. After consulting a key for the Junonia of French Guiana that relies on colour patterns and morphological traits [23], Borchers and Marcus [24] suggested that the two light-coloured Argentinian populations correspond to J. genoveva and either a genetically disparate population of the same species or an undescribed cryptic species while the dark-coloured population corresponds to J. wahlbergi. The results of this study suggest that there is no close affinity between the Argentinian forms and J. genoveva or J. wahlbergi from French Guiana, so the tentatively assigned scientific names should be revised. We refer to the 2 light-coloured forms as J. “hilaris” and J. sp. affin “hilaris” and the dark-coloured form as J. “flirtea.

The samples from French Guiana and the Caribbean are divided into two major clusters (Figure 5). The first cluster is composed of RAF Populations 1, 3A, 3B, 3C, 4, and 5. Overall, this cluster appears to be strongly associated with forms of Junonia whose larval host plants are in order Lamiales: J. litoralis from French Guiana and J. neildi from the Caribbean (both use larval host black mangrove, Avicennia germinans), J. genoveva from French Guiana (larval host Hyptis atrorubens), and J. zonalis from the Caribbean (larval hosts Stachytarpheta jamaicensis, S. urticifolia, and Lippia nodiflora) [23, 99]. Since all Junonia from the Caribbean feed on plants in the Lamiales, this cluster contains the bulk (19/30) of Caribbean specimens in our analysis. In contrast, the population cluster that includes Populations 2A and 2B contains the bulk of specimens from species whose larval host plants are in the order Scrophulariales: J. coenia, J. divaricata, J. evarete, and J. wahlbergi. Junonia divaricata and J. evarete use Utricularia hispida as their larval host, while J. wahlbergi and J. evarete use Agalinis hispidula [23]. Junonia coenia from North America feeds on a wide variety of larval hosts in the order Scrophulariales including several Agalinis species [18]. Sharing larval host plants may facilitate habitat overlap among extant forms and the cooccurrence of organisms is a necessary precondition for interspecific hybridization and gene flow to take place. The overall congruency between the population clustering based on RAF genotyping and host plant use supports the hypothesis of L. Brévignon and C. Brévignon [23] that host plant use defines two major lineages within New World Junonia. However, it should also be noted that there are exceptions to this congruence: 1 specimen of J. evarete clustered with the Lamiales feeders while 12 specimens of J. genoveva, 5 specimens of J. zonalis, and 5 specimens of J. neildi clustered with the Scrophulariales feeders (Figure 5). This might be interpreted as evidence for some gene flow between these two lineages.

Some species of Junonia were not readily distinguishable from one another in the Structure analysis of the RAF data. Junonia coenia, J. divaricata, J. evarete, and J. wahlbergi are all associated with the same RAF population cluster (2A-2B). We suspect that our inability to distinguish between these forms is due, at least in part, to artefact because these four species were represented by the smallest number of individuals in the RAF Structure analysis (between 2 and 10 individuals sampled depending on the species, Table 1). Thus, the statistical power of the algorithm is poor for these species [80]. With additional sampling of these forms a more robust analysis of population structure for these Junonia species will be possible. While we were not able to reliably separate J. zonalis and J. neildi from each other, we are able to conclude, based on the available data, that these Caribbean formsappear to be genetically differentiated from the taxa that occur in French Guiana with the vast majority of Caribbean specimens assigned to Populations 2B and 3C (Figure 5). This supports the taxonomic hypothesis of L. Brévignon and C. Brévignon [23], which elevated these taxa to full species. It also explains why it has been so challenging to apply taxonomic names based on South American types to forms found in the Caribbean [22, 2628].

The other pattern emerging from the analysis of RAF genotypes is for individuals of one species to be spread across multiple RAF populations (Figure 5). Junonia litoralis from French Guiana has several individuals assigned to the same RAF population (Population 1), but that population also includes individuals from at least two other species. Furthermore, other J. litoralis are assigned to a different population cluster (Populations 2A and 2B). In the most extreme case of distribution among RAF populations, J. genoveva from French Guiana has individuals assigned to 6 RAF populations (spread across two major population clusters with apparently little gene flow between the clusters). J. zonalis from the French Antilles is similarly distributed across 4 populations while J. neildi is assigned to at least 2 populations.

This is similar to what has been observed previously for the light buckeye, Junonia “hilaris” from Argentina [24] and has been replicated in the current analysis of RAF data, which divided this species into two separate populations (Populations 6A and 6B). Previously, we suggested that the presence of two very phenotypically similar, but genetically distinct, populations of Junonia existing simultaneously in Buenos Aires, Argentina, may be due to mass migrations of individuals from geographically disparate areas [24]. Mass migration of Junonia is a phenomenon that has been documented in Argentina and elsewhere in the New World [100, 101]. However, mass migration of Junonia has not been observed in French Guiana (C. Brévignon, pers. com.) and mass migration followed by hybridization between forms would tend to homogenize the genotypes of the interacting populations over time. While some of the named Junonia taxa from French Guiana and the French Antilles have been described relatively recently [23, 4749], all of these forms have been observed in the region for decades (C. Brévignon, pers. com.) and in some cases for centuries [25, 102].

Many different forms of New World Junonia tested in lab crosses are interfertile and produce viable fertile hybrids [3941], but many of these interfertile forms are separated geographically or by habitat preference (see below) and would have limited contact in the wild. Preliminary attempts at some interspecific pairings of sympatric Junonia from French Guiana have been unsuccessful (C. Brévignon, pers. com.). Wolbachia bacterial infections can prevent otherwise genetically compatible insects from producing viable offspring, blocking gene flow [103], perhaps contributing to the pattern we see in Junonia. Several species of Junonia have been tested for Wolbachia (including J. evarete from Panama), but thus far infections have only been detected in Asian J. almana [104, 105]. Junonia species have characteristic male genitalia [23, 24], but there is no evidence of male-female genitalic incompatibility as has been reported in some pairs of snail and moth sister species [106, 107]. Also, Junonia species have been observed to engage in courtship flights with heterospecifics in the wild [108] and wild-caught individuals that appear to be of hybrid origin have been identified based on their RAF genotypes [24] (Figure 5). This suggests that if mass migration were sufficiently common to bring individuals from geographically distinct populations into contact at a specific site, it would very quickly eliminate most of the genetic population structure in Junonia at that locality unless there is assortative mating between forms.

In addition to explaining patterns of genetic overlap among Junonia RAF populations, we have to explain the possible subdivision of some Junonia species across multiple populations (e.g., J. genoveva in French Guiana with 6 different RAF subpopulations, Figure 5). Frequently, individuals from one named taxon collected from a single locality were assigned to different RAF subpopulations. An alternative cause of the complex genetic population structure found in some forms of Junonia is that these species, which are defined on the basis of morphology and colour patterns, may include races that are specializing on different larval host plants. Host plant specialization is a widespread mechanism for population differentiation causing rapid evolution of adaptive traits for feeding on new hosts and for assortative mating to maintain favourable combinations of traits [109113]. In some cases, it has been suggested that this has been a driver for reproductive isolation and incipient speciation in many insects [113117].

Most New World Junonia are currently known to only feed on a single species of larval host plant in the wild [23, 49, 97]. However, J. coenia feeds on many alternative hosts [18]. Under artificial conditions, many varieties of Junonia larvae can be reared on alternate host plants or on artificial diets containing alternate host plant leaves ([40, 118] and Jeffrey M. Marcus, pers. observation). When presented with several alternative hosts, female J. coenia choose to oviposit on the same primary host plant used by the wild population from which the female was derived [21, 119]. If additional larval host plants for South American forms of Junonia exist, this may explain some of the extensive population structure seen in some species such as J. genoveva. New World Junonia host plants contain iridoid glycoside secondary compounds, the presence of which may be a necessary precondition for use as a host by Junonia [120]. This may help identify possible additional larval host plants for South American forms of Junonia.

A variety of mechanisms permit assortative mating and allow genetically distinct but reproductively compatible populations of species to persist in the same habitat. Habitat partitioning allows individuals from different populations to use different portions of available habitat, thereby making it less likely that they will interact and mate [121, 122]. In some North American habitats where multiple Junonia taxa occur, one of the authors (Jeffrey M. Marcus) has observed differences in habitat use. In mangrove swamps along the coast of Florida, USA, most J. coenia males appear to patrol mating territories in clearings without trees or other vertical habitat structure. In contrast, males of a form that resembles J. neildi (the forms have yet to be compared genetically) and whose larvae feed on black mangrove trees (Avicennia germinans) tolerate more vertical structure and establish mating territories in close proximity to their larval host plants. Similarly, in coastal dune habitats in south Texas, USA, J. coenia males establish mating territories on the foreshore between the sand dunes and the water line. A few meters away, males of the darkly pigmented J. “nigrosuffusa” (taxonomic affinities uncertain) appear to establish mating territories in the interdune and slack areas between the sand dunes. It is not clear whether these North American Junonia habitat preferences are due to preferences for abiotic conditions of the microhabitat itself or whether the presence or relative abundance of preferred larval host plants for each form in the favoured microhabitats is the driver of habitat preference [123]. In French Guiana, the presence of different Junonia taxa appears to be closely tied to the abundance and phenology of larval host plants [23]. Whether there are similar patterns of microhabitat subdivision among cooccurring Junonia species elsewhere is unknown.

A second mechanism for assortative mating is for different forms to become reproductively active at different times, reproducing in different years [124], at different times of the year [125], or at different times of the day [107, 126]. This reduces the likelihood of interspecific mating and permits the continued coexistence of allochronic species. There are no known diurnal differences in habitat usage among forms of Junonia, but there are seasonal differences that may contribute to an allochronic mechanism for persistence. In French Guiana, the foliage of Junonia larval host plants in the Lamiales is persistent while the foliage of larval hosts in the Scrophulariales deteriorates quickly during the dry season. The flight times of adults of J. divaricata, J. evarete, and J. wahlbergi coincide temporally with each other and with the presence of larval hosts in the Scrophulariales while the flight times of J. genoveva, J. litoralis, J. neildi, and J. zonalis are less seasonally restricted [23]. This difference in phenology may contribute to the continued distinctiveness of the two major Junonia lineages (clusters of RAF populations) in French Guiana and the French Antilles (Figure 5) but do not present an obvious mechanism for maintaining distinctive forms or species within a cluster.

A further possible mechanism may be differences in the amount or chemical composition of the pheromone or combination of pheromones used in the mating systems of different strains or species, which may allow individuals to establish a preference for other members of their own species [127]. Pheromones may differ because of intrinsic genetic differences between strains [128] or because of differences in the availability of pheromone precursors in the host plants used by different strains [129]. Unfortunately, nothing is currently known about Junonia pheromone use or composition. Other characteristics of mating systems may also contribute to the assortative mating between strains or species including vocalizations, displays of colour, and physical interactions between sexes [130, 131]. Of particular interest in this regard are several characteristics of the courtship flights that are known to differ among the North American and Caribbean forms of Junonia [26, 108]. Variation in courtship flight patterns in other Junonia forms has not yet been documented. There are also differences in the colour patterns of New World Junonia species [22, 23, 97], but any roles that these colour pattern differences play in the mating systems are also undocumented.

Operationally, we use the isolation species concept that defines species as systems of populations such that gene exchange between these systems is limited or prevented by one or more reproductive isolating mechanisms [45, 46].

5. Conclusions

While the molecular tools employed here cannot yet distinguish between all named forms of Junonia we are getting much closer to having a set of reliable molecular markers for defining groups of populations within which genetic exchange is extensive and between which genetic exchange is limited, providing a means by which we can begin to distinguish between species. While COI barcodes are of limited utility, nuclear wingless sequences and RAF genotyping are effective at identifying some individual species of Junonia and have been very helpful in examining the relationships of Junonia forms from different geographic regions. Using these tools, we have determined that, in spite of phenotypic similarities, Junonia from the French Antilles, French Guiana, and Argentina are genetically distinct from one another and that different species likely occur in each region. Junonia populations also appear to cluster according to larval host plant use, supporting the hypothesis that there are two Junonia lineages: one which feeds primarily on plants in the order Scrophulariales and the other which uses larval host plants in the order Lamiales. The rapid growth in our knowledge of the natural and evolutionary history of New World Junonia in combination with the powerful experimental tools that are available for use in these organisms shows much promise in making this group an excellent model for the study of processes of speciation, host plant adaptation, and the evolution and development of colour pattern phenotypes.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors’ sincere thanks are due to Christian Brévignon for providing the Junonia samples analyzed in this paper and for his extensive and insightful comments on an earlier draft of this paper. The authors also benefitted greatly from the helpful comments of Martin Villet, Niklas Wahlberg, and two anonymous reviewers. Thanks are due to Roohollah Abbasi, Kahlia Beaudette, Ashley Haverstick, Moiz Kapasi, Bonnie McCullagh, Jacob Miller, Melissa Peters, and Stephanie Rozbacher for assistance in the laboratory and for helpful comments during the preparation of this paper. They would like to thank Margaret Docker for the use of the gel imager. Thanks are due to Rodolphe Rougerie, Evgeny Zakharov, and Suresh Naik for their assistance in obtaining access to the Junonia DNA extractions processed at the Canadian Centre for DNA Barcoding at the University of Guelph. Finally, they thank Naomi Rowland from the Biotechnology Core Facility of Western Kentucky University and the staff of the University Core DNA Services Lab of the University of Calgary for processing of RAF and sequencing samples, respectively. Support for this project was provided by a University of Manitoba Undergraduate Research Award and an NSERC-USRA (to Amber P. Gemmell), in addition to an NSERC Discovery Grant RGPIN386337-2011 and a Canada Foundation for Innovation Award (to Jeffrey M. Marcus).