# Molecular Population Structure ofJunonia Butterflies from French Guiana, Guadeloupe, and Martinique

Authors: Amber P. Gemmell; Tanja E. Borchers; Jeffrey M. Marcus
Journal: Psyche
 (2014)

---

## Abstract

Up to 9 described species ofJunonia 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.


---

## Body

## 1. Introduction

Buckeye butterflies, genusJunonia (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 [2–9]. Experimental tools to manipulate gene expression developed inJunoniaare broadly applicable across the Lepidoptera [10–14].Junonia has also been used in studies of insect endocrinology [15–17] and has been an important system for examining the evolution of larval host plant preference and tolerance to host plant toxins [18–21].Junonia butterflies are found throughout the Old and New World tropics. In the Western Hemisphere, forms ofJunonia occur from southern Canada to Tierra del Fuego [22–24] and have a complicated taxonomic history. In 1775, Cramer [25] identified and described two similar species ofJunonia,J. evarete andJ. 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 toJ. evarete [29–33]. This is a result of the geographical [29, 34] and seasonal [35] variability ofJunonia and the fact that someJunonia forms closely resemble one another [36, 37]. In addition, differentJunonia forms share identical karyotypes (N
=
31) [26, 38] and are capable of hybridization and the production of fertile offspring [39–41], further complicating the process of assigning names toJunonia specimens. These features make applying the biological species concept [42], phylogenetic species concept [43], or morphospecies concept [44] very difficult inJunonia. 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 inJunonia will be of great importance in clarifyingJunonia taxonomy.Authorities who favoured the two-species hypothesis inJunonia called the larger formJ. genoveva and the smaller formJ. evarete. In 1985, Turner and Parnell [26], using specimens from Jamaica and Florida, USA, verified the existence of twoJunonia 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 nowJ. evarete and the smaller species wasJ. 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 andJunonia specimens from many localities in South America. Using this reference material, Neild [22] reversed Turner and Parnell [26] so that the larger species was againJ. genoveva and the smaller species wasJ. evarete. Neild [22] also designated new types forJ. evarete andJ. genoveva to facilitate future taxonomic work.Recently, L. Brévignon and C. Brévignon [23, 47, 48] identified 5Junonia species (J. evarete,J. genoveva,J. wahlbergi,J. litoralis, andJ. divaricata) from French Guiana. This represents the most diverse assemblage ofJunonia in the New World. There are two forms ofJunoniaknown from the Caribbean Islands, “zonalis” and “neildi,” which were initially recognized as subspecies of mainlandJ. evarete andJ. genoveva, respectively [49] and later as two distinct species:J. zonalis andJ. neildi [23]. In recognizingJ. neildi andJ. zonalis as distinct species, L. Brévignon and C. Brévignon [23] restricted the use of the species epithetsJ. evarete andJ. genovevato mainland Central and South American forms. If it were confirmed that the Caribbean forms are actually distinct species with respect toJunonia 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 additionalJunonia species,J. coenia from North America andJ. vestina from the Andes mountains of South America, for a current total of up to 9 species of New WorldJunonia.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 genusJunonia [50, 51] (some authorities had previously placed these species in the related genusPrecis) [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 WorldJunonia [23, 54, 55], which have better taxon sampling, have focused entirely on the mitochondrialcytochrome 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 WorldJunonia 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 inJunonia is unknown because only mitochondrial markers were considered. What is apparent is that there are two very divergentCOI haplotype groups (4% sequence divergence between them) present in New WorldJunonia: 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 WorldJunonia taxa using molecular markers employed a combination of mitochondrial and nuclear markers to examine populations ofJunonia in Buenos Aires, Argentina. Borchers and Marcus [24] used DNA sequences from the nuclearwingless gene and anonymous nuclear loci identified by Randomly Amplified Fingerprinting (RAF) (a technique used to assess genetic diversity within populations [58–60] and gene flow between populations [61]) in addition to sequences from the mitochondrialCOI gene. They identified 3 distinct populations ofJunonia from Buenos Aires: one population with dark-coloured wings referred to asJ. evarete flirtea [62] which Borchers and Marcus [24] suggested may correspond toJ. wahlbergi and 2 light-coloured populations that correspond toJ. genoveva hilaris and either a genetically disparate population ofJ. genoveva hilaris or an undescribed crypticJunonia species. However, the relationship of these Argentinian forms withJ. evarete andJ. genovevafrom Suriname and French Guiana is not known, so we refer to them asJ. “flirtea” andJ. “hilaris.”In the current study, we extend the genetic tools that were employed by Borchers and Marcus [24] toJunonia 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 studiedJunonia from Argentina [24].


## 2. Materials and Methods

### 2.1. Specimens and DNA Preparation

A total of 104Junonia 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.Table 1
Number ofJunoniaspecimens included in this study either entirely processed in our laboratory or extracted at the University of Guelph and sent to us for further study.
Species and locality
DNA extracted in our laboratory
DNA extracted by Guelph and whole genome amplified by our laboratory
J. coenia,Florida, USA
0
2
J. divaricata,French Guiana
0
5
J. evarete,French Guiana
0
6
J. genoveva, French Guiana
32
6
J. litoralis, French Guiana
8
4
J. neildi, Guadeloupe
6
2
J. neildi, Martinique
2
3
J. wahlbergi, French Guiana
0
10
J. zonalis, Guadeloupe
7
2
J. zonalis, Martinique
7
2Only 10μL aliquots of DNA were available for the 42Junonia 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. MitochondrialCytochrome 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).Table 2
Primer sequences used incytochrome oxidase I (COI) and wingless(wg) PCR reactions.
Primer name
Sequence
COI
Nancy
5′CCCGGTAAAATTAAAATATAAACTTC3′
LCO1490
5′GGTCAACAAATCATAAAGATATTGG3′
HCO2198
5′TAAACTTCAGGGTGACC AAAAAATCA3′
M13-uniminibarF1
5′GTAAAACGACGGCCAGTGGAAAATCATAATGAAGGCATGAGC3′
M13-uniminibarR1
5′GGAAACAGCTATGACCATGTCCACTAATCACAARGATATTGGTAC3′
miniCOIF2
5′ATACTATTGTTACAGCCTCATGC3′
miniCOIR2
5′TGTTGTAATAAAATTAATAGCTCC3′
miniCOIF3
5′CCCCACTTTCATCTAATATTGC3′
wg
lepwg1
5′GARTGYAARTGYCAYGGYATGTCTGG3′
lepwg2
5′ACTNCGCRCACCATGGAATGTRCA3′
miniwgF
5′ATCGCGGGTCATGATGCCTAATACG3′
miniwgR
5′GTTCTTTTCGCAGAAACCCGGTGAAC3′
T7-miniwgR
5′TAATACGACTCACTATAGGGGTTCTTTTCGCAGAAACCCGGTGAAC3′Samples that failed to amplify with LCO1490 and HCO2198 were reamplified with M13-uniminibarF1 (miniCOIF) and M13-uniminibarR1 (miniCOIR) (Table2) [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 theJunonia 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. NuclearWingless 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 inJunonia [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 theJunonia 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.Table 3
Expected sequence length of trimmed PCR products (primers removed) for each primer pair.
Primer pair
Trimmed sequence length (base pairs)
LCO1490/Nancy
725
LCO1490/HCO2198
658
M13-uniminibarF1/M13-uniminibarR1
153
mCOIF2/mCOIR2
292
mCOIF2/HCO2198
520
mCOIF3/HCO2198
295
Lepwg1/Lepwg2
402
miniwgF/miniwgR
137


### 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. MitochondrialCytochrome 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 resequencedCOI 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 publishedJunonia COI sequences were included in this phylogenetic analysis [23, 24, 51, 57, 70–74]. We also conducted a maximum likelihood bootstrap analysis of this dataset (500 fast addition replicates, collapsing all notes with frequency less than 50%). The alignedCOI FASTA sequences generated by this study along with 22 previously published ArgentinianJunonia 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. NuclearWingless Analysis

For theJunonia species sequenced in this study and ArgentinianJunonia 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 probablewingless 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 (ExactG 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 thewingless data since, unlike GENEPOP [78], Structure does not require thea priori assignment of individuals to specific subpopulations. Population structure exhibited bywingless 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 88wingless sequences. The maximum log likelihood (ln
⁡
⁡
P
(
D
)) for the 10 replicate searches for each population model was used to calculate the posterior probability (P
(
K
=
n
)) of each population model. Haplotype networks ofwingless alleles were constructed in the same manner asCOI except that PHASE output identifying the most likelywingless 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 ArgentinianJunonia [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 thewingless 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 K = 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 forCOI andwingless.


## 2.1. Specimens and DNA Preparation

A total of 104Junonia 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.Table 1
Number ofJunoniaspecimens included in this study either entirely processed in our laboratory or extracted at the University of Guelph and sent to us for further study.
Species and locality
DNA extracted in our laboratory
DNA extracted by Guelph and whole genome amplified by our laboratory
J. coenia,Florida, USA
0
2
J. divaricata,French Guiana
0
5
J. evarete,French Guiana
0
6
J. genoveva, French Guiana
32
6
J. litoralis, French Guiana
8
4
J. neildi, Guadeloupe
6
2
J. neildi, Martinique
2
3
J. wahlbergi, French Guiana
0
10
J. zonalis, Guadeloupe
7
2
J. zonalis, Martinique
7
2Only 10μL aliquots of DNA were available for the 42Junonia 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. MitochondrialCytochrome 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).Table 2
Primer sequences used incytochrome oxidase I (COI) and wingless(wg) PCR reactions.
Primer name
Sequence
COI
Nancy
5′CCCGGTAAAATTAAAATATAAACTTC3′
LCO1490
5′GGTCAACAAATCATAAAGATATTGG3′
HCO2198
5′TAAACTTCAGGGTGACC AAAAAATCA3′
M13-uniminibarF1
5′GTAAAACGACGGCCAGTGGAAAATCATAATGAAGGCATGAGC3′
M13-uniminibarR1
5′GGAAACAGCTATGACCATGTCCACTAATCACAARGATATTGGTAC3′
miniCOIF2
5′ATACTATTGTTACAGCCTCATGC3′
miniCOIR2
5′TGTTGTAATAAAATTAATAGCTCC3′
miniCOIF3
5′CCCCACTTTCATCTAATATTGC3′
wg
lepwg1
5′GARTGYAARTGYCAYGGYATGTCTGG3′
lepwg2
5′ACTNCGCRCACCATGGAATGTRCA3′
miniwgF
5′ATCGCGGGTCATGATGCCTAATACG3′
miniwgR
5′GTTCTTTTCGCAGAAACCCGGTGAAC3′
T7-miniwgR
5′TAATACGACTCACTATAGGGGTTCTTTTCGCAGAAACCCGGTGAAC3′Samples that failed to amplify with LCO1490 and HCO2198 were reamplified with M13-uniminibarF1 (miniCOIF) and M13-uniminibarR1 (miniCOIR) (Table2) [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 theJunonia 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. NuclearWingless 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 inJunonia [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 theJunonia 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.Table 3
Expected sequence length of trimmed PCR products (primers removed) for each primer pair.
Primer pair
Trimmed sequence length (base pairs)
LCO1490/Nancy
725
LCO1490/HCO2198
658
M13-uniminibarF1/M13-uniminibarR1
153
mCOIF2/mCOIR2
292
mCOIF2/HCO2198
520
mCOIF3/HCO2198
295
Lepwg1/Lepwg2
402
miniwgF/miniwgR
137


## 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. MitochondrialCytochrome 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 resequencedCOI 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 publishedJunonia COI sequences were included in this phylogenetic analysis [23, 24, 51, 57, 70–74]. We also conducted a maximum likelihood bootstrap analysis of this dataset (500 fast addition replicates, collapsing all notes with frequency less than 50%). The alignedCOI FASTA sequences generated by this study along with 22 previously published ArgentinianJunonia 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. NuclearWingless Analysis

For theJunonia species sequenced in this study and ArgentinianJunonia 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 probablewingless 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 (ExactG 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 thewingless data since, unlike GENEPOP [78], Structure does not require thea priori assignment of individuals to specific subpopulations. Population structure exhibited bywingless 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 88wingless sequences. The maximum log likelihood (ln
⁡
⁡
P
(
D
)) for the 10 replicate searches for each population model was used to calculate the posterior probability (P
(
K
=
n
)) of each population model. Haplotype networks ofwingless alleles were constructed in the same manner asCOI except that PHASE output identifying the most likelywingless 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 ArgentinianJunonia [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 thewingless 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 K = 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 forCOI andwingless.


## 3. Results

### 3.1. MitochondrialCytochrome Oxidase I Results

NewCOI DNA sequences generated by this project were deposited in Genbank (67 accessions, numbers KJ469059–KJ469126), with the exception of specimens with onlyCOI minibarcode sequence fragments [65], which were submitted to the DNA Databank of Japan (DDBJ, 5 accessions AB935341–AB935345). FullCOI 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, someCOI sequence was recovered from 90 of the 104 specimens.Analysis of theCOI sequences produced a maximum likelihood phylogenetic tree (Figure 1). As previously reported [24, 54], there are two distinct mitochondrial haplotype groups in New WorldJunonia. 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 ofJunonia appear to be associated with only one haplotype group (group A: the South American formsJ. “flirtea” andJ. vestina; group B: the North American formsJ. coeniaandJ. “nigrosuffusa”). All otherJunonia species (J. divaricata,J. evarete,J. genoveva,J. “hilaris,” J. litoralis,J. neildi,J. wahlbergi, andJ. 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.Figure 1
Phylogenetic tree depicting the two distinct mitochondrialCOI haplotype groups found among New WorldJunonia. Group A haplotypes occur primarily in South America and the Caribbean. Group B haplotypes occur throughout the Western Hemisphere but are most common in North and Central America. MostJunonia species include individuals with both haplotypes. OnlyJ. “flirtea” andJ. vestina possess exclusively group A haplotypes, while onlyJ. coenia andJ. evarete nigrosuffusa possess exclusively group B haplotypes. In this figure, the entire haplotype Group B clade has been translated horizontally so that the tree fits better on the page. The dotted portion of the line connecting this clade to the tree is not included in calculations of sequence distances.Haplotype networks generated using complete barcode fragment mitochondrialCOI haplotypes. Circles are scaled to represent the number of individuals that contain a specificCOI haplotype. Divisions and colours of circles in (a) reflect geography and species associated with eachCOI haplotype. French GuiananJ. genoveva is orange, French GuiananJ. litoralis is yellow, French GuiananJ. wahlbergi is red, French GuiananJ. evarete is pale orange, French GuiananJ. divaricata is pale yellow, AmericanJ. coenia is grey, MartiniquanJ. zonalis is purple, MartiniquanJ. neildi is pink, GuadeloupeanJ. zonalis is dark blue, GuadeloupeanJ. neildi is light blue, ArgentinianJ. “hilaris” is light green, and ArgentinianJ. “flirtea” is dark green. Colours of circles in (b) reflect the RAF population associated with eachCOI haplotype. Population 1 is yellow, Population 2A is dark blue, Population 2B is light blue, Population 3A is dark purple, Population 3B is light purple, Population 3C is pink, Population 4 is red, Population 5 is orange, Population 6A is light green, Population 6B is medium green, and Population 6C is dark green. All 11 RAF populations are represented in group A haplotypes. Group B haplotypes are comprised of RAF Populations 2A, 2B, 3C, 5, and 6B.
(a)
(b)Haplotype networks constructed from those samples with partial mCOIF2/HCO fragmentCOI sequences. The network in (a) reflects geography and species while the network in (b) reflects RAF population assignment. The key to the colours is the same as in Figure 2.
(a)
(b)Junonia litoralis from French Guiana andJ. neildi from the Caribbean, both of which feed on black mangrove (Avicennia germinans) as larvae, include individuals with mitochondria from bothCOI haplotype groups.Junonia genoveva from French Guiana, which has often been considered conspecific withJ. neildi from the Caribbean, also includes individuals that carry haplotypes in groups A and B.Junonia evarete from French Guiana have exclusively group ACOI haplotypes, butJ. zonalis from the Caribbean,which has sometimes been considered conspecific withJ. 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 fromCOI 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 ofJ. 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 theCOI 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 allJ. wahlbergi and almost allJ. evarete from French Guiana share a single group A genotype that is rare in otherJunoniaspecies. Similarly, the vast majority ofJ. genoveva from French Guiana possess one of 3 disparate genotypes withinCOI haplotype group A (Figures 2(a) and 3(a)) that are rare in all otherJunonia species. Most of the remainingJ. 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 ofJunonia there are no abundantCOI genotypes that are diagnostic for particular species.


### 3.2. NuclearWingless Results

Short 137 bpwingless 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 22wingless sequences from ArgentinianJunonia specimens [24]. Following CLUSTALW alignments of the sequences of the mini-wingless PCR products with the Argentinianwingless 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 35wingless alleles. The most probable allelic combinations for eachJunonia 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. SeparatingJunonia populations solely by geography (P
=
0.00017) or species (P
=
0.0000) and by both geography and species (P
=
0.0000) suggests significant genetic differentiation and distinctwingless allele distributions associated with these two factors. However, separating individualJunonia specimens using mitochondrialCOI haplotype alone to define populations shows no statistically significant distinct distribution ofwingless alleles associated with mitochondrial haplotypes (P
=
0.776). However, when specimens were categorized by geography andCOI haplotype (P
=
0.00014); species andCOI haplotype (P
=
0.0000); or geography, species, andCOI haplotype (P
=
0.0000), thewingless 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 ofwingless alleles.Thewingless 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 (P
(
K
=
n
)) based on thewingless allelic data is K
=
1 (all samples belonging to a single population). Haplotype networks ofwingless alleles show that one allele (a9) that is common inJ. genoveva populations from French Guiana (50%) is much rarer inJunonia from both Martinique and Guadeloupe (15%) (Figure 4). A secondJunonia wingless allele (a7) is fairly common in both French Guianan (38%) and French Antillean (54%) populations.Table 4
InferringK, the number of populations, testing for 1–15 subpopulations in STRUCTURE, for Argentinian, French Guianan, Guadeloupean, and Martiniquan Junonia wingless data.
K (number of populations in model)
Medianln
⁡
[
P
(
D
)
]
PR (K
=
n) (posterior probability of the model)
1
−42.70
0.972
2
−46.54
0.021
3
−56.03
1.587
E
-
06
4
−53.68
1.664
E
-
05
5
−136.83
1.287
E
-
41
6
−135.59
4.426
E
-
41
7
−173.16
2.136
E
-
57
8
−130.03
1.156
E
-
38
9
−149.33
4.797
E
-
47
10
−61.83
4.806
E
-
09
11
−59.65
4.251
E
-
08
12
−55.59
2.465
E
-
06
13
−53.93
1.296
E
-
05
14
−49.28
0.001
15
−47.82
0.006Haplotype networks generated usingwingless alleles. Circles are scaled to represent the number of individuals that contain a specificwingless allele with the exception of alleles a7 and a9 (if scaled proportionately, would be 50 and 61 times larger, resp., than shown). Colours for both (a) and (b) are as described in Figure 2. (a) Divisions and colours of circles reflect geography and species associated with eachwingless allele. Allele a9 is much rarer in the Caribbean populations than it is in the mainland populations. 50% of a9 allele is comprised ofJ. genoveva. The majority ofwingless alleles found in ArgentinianJunonia are allele a9 or its derivatives. (b) Divisions and colours of circles reflect the RAF populations associated with eachwingless allele. All 11 RAF populations carry the most common allele a9. Individuals from 10 of the 11 RAF populations carry the other common allele, a7.
(a)
(b)


### 3.3. Randomly Amplified Fingerprinting Results

RAF fragments for the RP2, RP4, and RP6 primers were recovered from all 104Junonia specimens in this study and 22 ArgentinianJunonia 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 CaribbeanJunonia butterflies. The results of this analysis (Table 6) showed that the model with the highest posterior probability (P
(
K
=
n
)
=
1) is K
=
8 (samples belonging to 8 separate populations). Curiously, when this analysis was repeated with the inclusion ofJunonia samples from Argentina, the model with the highest posterior probability (P
(
K
=
n
)
=
1) was K
=
6. However, the Structure analysis did not distribute the Argentinian samples among the populations previously established for French Guiana and the Caribbean. Instead, all of theJunonia 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 ArgentinianJunonia, 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 ArgentinianJunonia [24]. When the Argentinian data set is analyzed in isolation both sets of. RAF loci show the model with the highest posterior probability (P
(
K
=
n
)
=
1) is K
=
3 (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).Table 5
Fragment sizes, in base pairs, of RAF amplification products for RP2, RP4, and RP6 primer.
RP2 (bp)
RP4 (bp)
RP6 (bp)
35
37
35
39
42
38
41
46
47
43
50
50
47
56
61
49
58
63
51
62
72
54
65
73
56
68
77
58
71
82
60
79
85
64
82
89
66
85
93
68
87
147
72
89
149
75
93
190
76
95
193
79
99
197
80
100
84
103
86
112
87
141
88
145
90
149
92
152
95
155
98
158
100
161
107
164
135
168
136
175
148
179
150
196
158
201
160
202
162
209
182
223
200
228
228
231
233
234
234
237
240
240
241
247
249
258
261
277
280
287
289
291
342
345
363
369
371
389
391
394
440
483Table 6
InferringK, the number of populations, testing for 1–15 subpopulations in STRUCTURE, for the Argentinian, French Guianan, Guadeloupean, and Martiniquan JunoniaRAF data.
K(number of populations in model)
Medianln
⁡
[
P
(
D
)
]
PR (K
=
n) (posterior probability of the model)
Without ArgentinianJunonia
1
−235.98
9.598
E
-
47
2
−178.72
7.078
E
-
22
3
−147.20
3.458
E
-
08
4
−143.08
2.129
E
-
06
5
−136.83
0.001
6
−135.59
0.004
7
−173.16
1.839
E
-
19
8
−130.03
0.995
9
−149.33
4.130
E
-
09
10
−152.54
1.659
E
-
10
11
−146.13
1.008
E
-
07
12
−160.62
5.136
E
-
14
13
−150.26
1.630
E
-
09
14
−157.84
8.279
E
-
13
15
−157.205
1.562
E
-
12
With ArgentinianJunonia
1
−384.19
3.286
E
-
78
2
−300.45
7.665
E
-
42
3
−237.34
1.963
E
-
14
4
−208.28
0.082
5
−207.97
0.112
6
−206.27
0.611
7
−214.73
0.000
8
−208.28
0.082
9
−299.36
2.291
E
-
41
10
−313.86
1.156
E
-
47
11
−342.71
3.398
E
-
60
12
−258.91
8.416
E
-
24
13
−207.97
0.112
14
−243.61
3.714
E
-
17
15
−270.82
5.658
E
-
29Summary of population genetic analysis of RAF nuclear genotyping. (a) Venn diagram illustrating gene flow occurring between populations. The 11Junonia populations are distributed into 3 clusters which show extensive genetic exchange within a cluster, but little or no genetic exchange between clusters. (b) CONTML analysis of the RAF allele frequency data for populations identified in Structure. Clusters identified in (a) tend to be adjacent to one another on the tree, which shows a very strong influence on larval host plant use. Most of the individuals in RAF Populations 2A, 2B, 6A, 6B, and 6C are from species that utilize larval host plants in the order Scrophulariales. Virtually all of theJunonia in the remaining RAF populations are from species that use larval host plants in the order Lamiales.
(a)
(b)Overall, there are 11 populations established by Structure analysis of RAF alleles (Figure5(a)). RAF Population 1 (yellow) primarily includes specimens ofJ. litoralis, but 2J. zonalis, 1J. genoveva, and 1J. evarete specimens also show genetic influence from this population.Population 2A (dark blue) includes all specimens ofJ. wahlbergi andJ. divaricata and all but one of theJ. evarete specimens. Population 2A also includes individuals from all of the other French Guianan and French AntilleanJunonia species. Most of the specimens in Population 2B (light blue) also show genetic influence from Population 2A, but there are 2 CaribbeanJ. zonalis specimens whose primary genetic influence is population 2B. Population 3A (dark purple) is comprised ofJ. genovevaandJ. zonalis specimens. Populations 3B (light purple) and 3C (pink) both includeJ. genoveva,J. neildi, andJ. zonalisspecimens. Population 4 (red) contains onlyJ. genoveva specimens. Population 5 (orange) is primarily composed ofJ. genovevaspecimens, although there is 1J. zonalis specimen that also shows influence from this population. Populations 6A-C are the same as the 3 ArgentinianJunonia 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 (Figure5(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 ofJ. zonalis andJ. neildi from the Caribbean andJ. genovevaandJ. litoralis from French Guiana comprise the remainder of this cluster (along with 1 specimen ofJ. 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 fromJ. evarete,J. wahlbergi, andJ. divaricata from French Guiana, although it also includes specimens from the other 2 French Guianan species, specimens ofJ. neildiandJ. zonalis from the Caribbean, andJ.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 (Figure5(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 withCOI andwingless allele frequency data produces a tree of identical topology with minuscule changes in branch length (not shown).


## 3.1. MitochondrialCytochrome Oxidase I Results

NewCOI DNA sequences generated by this project were deposited in Genbank (67 accessions, numbers KJ469059–KJ469126), with the exception of specimens with onlyCOI minibarcode sequence fragments [65], which were submitted to the DNA Databank of Japan (DDBJ, 5 accessions AB935341–AB935345). FullCOI 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, someCOI sequence was recovered from 90 of the 104 specimens.Analysis of theCOI sequences produced a maximum likelihood phylogenetic tree (Figure 1). As previously reported [24, 54], there are two distinct mitochondrial haplotype groups in New WorldJunonia. 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 ofJunonia appear to be associated with only one haplotype group (group A: the South American formsJ. “flirtea” andJ. vestina; group B: the North American formsJ. coeniaandJ. “nigrosuffusa”). All otherJunonia species (J. divaricata,J. evarete,J. genoveva,J. “hilaris,” J. litoralis,J. neildi,J. wahlbergi, andJ. 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.Figure 1
Phylogenetic tree depicting the two distinct mitochondrialCOI haplotype groups found among New WorldJunonia. Group A haplotypes occur primarily in South America and the Caribbean. Group B haplotypes occur throughout the Western Hemisphere but are most common in North and Central America. MostJunonia species include individuals with both haplotypes. OnlyJ. “flirtea” andJ. vestina possess exclusively group A haplotypes, while onlyJ. coenia andJ. evarete nigrosuffusa possess exclusively group B haplotypes. In this figure, the entire haplotype Group B clade has been translated horizontally so that the tree fits better on the page. The dotted portion of the line connecting this clade to the tree is not included in calculations of sequence distances.Haplotype networks generated using complete barcode fragment mitochondrialCOI haplotypes. Circles are scaled to represent the number of individuals that contain a specificCOI haplotype. Divisions and colours of circles in (a) reflect geography and species associated with eachCOI haplotype. French GuiananJ. genoveva is orange, French GuiananJ. litoralis is yellow, French GuiananJ. wahlbergi is red, French GuiananJ. evarete is pale orange, French GuiananJ. divaricata is pale yellow, AmericanJ. coenia is grey, MartiniquanJ. zonalis is purple, MartiniquanJ. neildi is pink, GuadeloupeanJ. zonalis is dark blue, GuadeloupeanJ. neildi is light blue, ArgentinianJ. “hilaris” is light green, and ArgentinianJ. “flirtea” is dark green. Colours of circles in (b) reflect the RAF population associated with eachCOI haplotype. Population 1 is yellow, Population 2A is dark blue, Population 2B is light blue, Population 3A is dark purple, Population 3B is light purple, Population 3C is pink, Population 4 is red, Population 5 is orange, Population 6A is light green, Population 6B is medium green, and Population 6C is dark green. All 11 RAF populations are represented in group A haplotypes. Group B haplotypes are comprised of RAF Populations 2A, 2B, 3C, 5, and 6B.
(a)
(b)Haplotype networks constructed from those samples with partial mCOIF2/HCO fragmentCOI sequences. The network in (a) reflects geography and species while the network in (b) reflects RAF population assignment. The key to the colours is the same as in Figure 2.
(a)
(b)Junonia litoralis from French Guiana andJ. neildi from the Caribbean, both of which feed on black mangrove (Avicennia germinans) as larvae, include individuals with mitochondria from bothCOI haplotype groups.Junonia genoveva from French Guiana, which has often been considered conspecific withJ. neildi from the Caribbean, also includes individuals that carry haplotypes in groups A and B.Junonia evarete from French Guiana have exclusively group ACOI haplotypes, butJ. zonalis from the Caribbean,which has sometimes been considered conspecific withJ. 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 fromCOI 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 ofJ. 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 theCOI 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 allJ. wahlbergi and almost allJ. evarete from French Guiana share a single group A genotype that is rare in otherJunoniaspecies. Similarly, the vast majority ofJ. genoveva from French Guiana possess one of 3 disparate genotypes withinCOI haplotype group A (Figures 2(a) and 3(a)) that are rare in all otherJunonia species. Most of the remainingJ. 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 ofJunonia there are no abundantCOI genotypes that are diagnostic for particular species.


## 3.2. NuclearWingless Results

Short 137 bpwingless 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 22wingless sequences from ArgentinianJunonia specimens [24]. Following CLUSTALW alignments of the sequences of the mini-wingless PCR products with the Argentinianwingless 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 35wingless alleles. The most probable allelic combinations for eachJunonia 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. SeparatingJunonia populations solely by geography (P
=
0.00017) or species (P
=
0.0000) and by both geography and species (P
=
0.0000) suggests significant genetic differentiation and distinctwingless allele distributions associated with these two factors. However, separating individualJunonia specimens using mitochondrialCOI haplotype alone to define populations shows no statistically significant distinct distribution ofwingless alleles associated with mitochondrial haplotypes (P
=
0.776). However, when specimens were categorized by geography andCOI haplotype (P
=
0.00014); species andCOI haplotype (P
=
0.0000); or geography, species, andCOI haplotype (P
=
0.0000), thewingless 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 ofwingless alleles.Thewingless 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 (P
(
K
=
n
)) based on thewingless allelic data is K
=
1 (all samples belonging to a single population). Haplotype networks ofwingless alleles show that one allele (a9) that is common inJ. genoveva populations from French Guiana (50%) is much rarer inJunonia from both Martinique and Guadeloupe (15%) (Figure 4). A secondJunonia wingless allele (a7) is fairly common in both French Guianan (38%) and French Antillean (54%) populations.Table 4
InferringK, the number of populations, testing for 1–15 subpopulations in STRUCTURE, for Argentinian, French Guianan, Guadeloupean, and Martiniquan Junonia wingless data.
K (number of populations in model)
Medianln
⁡
[
P
(
D
)
]
PR (K
=
n) (posterior probability of the model)
1
−42.70
0.972
2
−46.54
0.021
3
−56.03
1.587
E
-
06
4
−53.68
1.664
E
-
05
5
−136.83
1.287
E
-
41
6
−135.59
4.426
E
-
41
7
−173.16
2.136
E
-
57
8
−130.03
1.156
E
-
38
9
−149.33
4.797
E
-
47
10
−61.83
4.806
E
-
09
11
−59.65
4.251
E
-
08
12
−55.59
2.465
E
-
06
13
−53.93
1.296
E
-
05
14
−49.28
0.001
15
−47.82
0.006Haplotype networks generated usingwingless alleles. Circles are scaled to represent the number of individuals that contain a specificwingless allele with the exception of alleles a7 and a9 (if scaled proportionately, would be 50 and 61 times larger, resp., than shown). Colours for both (a) and (b) are as described in Figure 2. (a) Divisions and colours of circles reflect geography and species associated with eachwingless allele. Allele a9 is much rarer in the Caribbean populations than it is in the mainland populations. 50% of a9 allele is comprised ofJ. genoveva. The majority ofwingless alleles found in ArgentinianJunonia are allele a9 or its derivatives. (b) Divisions and colours of circles reflect the RAF populations associated with eachwingless allele. All 11 RAF populations carry the most common allele a9. Individuals from 10 of the 11 RAF populations carry the other common allele, a7.
(a)
(b)


## 3.3. Randomly Amplified Fingerprinting Results

RAF fragments for the RP2, RP4, and RP6 primers were recovered from all 104Junonia specimens in this study and 22 ArgentinianJunonia 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 CaribbeanJunonia butterflies. The results of this analysis (Table 6) showed that the model with the highest posterior probability (P
(
K
=
n
)
=
1) is K
=
8 (samples belonging to 8 separate populations). Curiously, when this analysis was repeated with the inclusion ofJunonia samples from Argentina, the model with the highest posterior probability (P
(
K
=
n
)
=
1) was K
=
6. However, the Structure analysis did not distribute the Argentinian samples among the populations previously established for French Guiana and the Caribbean. Instead, all of theJunonia 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 ArgentinianJunonia, 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 ArgentinianJunonia [24]. When the Argentinian data set is analyzed in isolation both sets of. RAF loci show the model with the highest posterior probability (P
(
K
=
n
)
=
1) is K
=
3 (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).Table 5
Fragment sizes, in base pairs, of RAF amplification products for RP2, RP4, and RP6 primer.
RP2 (bp)
RP4 (bp)
RP6 (bp)
35
37
35
39
42
38
41
46
47
43
50
50
47
56
61
49
58
63
51
62
72
54
65
73
56
68
77
58
71
82
60
79
85
64
82
89
66
85
93
68
87
147
72
89
149
75
93
190
76
95
193
79
99
197
80
100
84
103
86
112
87
141
88
145
90
149
92
152
95
155
98
158
100
161
107
164
135
168
136
175
148
179
150
196
158
201
160
202
162
209
182
223
200
228
228
231
233
234
234
237
240
240
241
247
249
258
261
277
280
287
289
291
342
345
363
369
371
389
391
394
440
483Table 6
InferringK, the number of populations, testing for 1–15 subpopulations in STRUCTURE, for the Argentinian, French Guianan, Guadeloupean, and Martiniquan JunoniaRAF data.
K(number of populations in model)
Medianln
⁡
[
P
(
D
)
]
PR (K
=
n) (posterior probability of the model)
Without ArgentinianJunonia
1
−235.98
9.598
E
-
47
2
−178.72
7.078
E
-
22
3
−147.20
3.458
E
-
08
4
−143.08
2.129
E
-
06
5
−136.83
0.001
6
−135.59
0.004
7
−173.16
1.839
E
-
19
8
−130.03
0.995
9
−149.33
4.130
E
-
09
10
−152.54
1.659
E
-
10
11
−146.13
1.008
E
-
07
12
−160.62
5.136
E
-
14
13
−150.26
1.630
E
-
09
14
−157.84
8.279
E
-
13
15
−157.205
1.562
E
-
12
With ArgentinianJunonia
1
−384.19
3.286
E
-
78
2
−300.45
7.665
E
-
42
3
−237.34
1.963
E
-
14
4
−208.28
0.082
5
−207.97
0.112
6
−206.27
0.611
7
−214.73
0.000
8
−208.28
0.082
9
−299.36
2.291
E
-
41
10
−313.86
1.156
E
-
47
11
−342.71
3.398
E
-
60
12
−258.91
8.416
E
-
24
13
−207.97
0.112
14
−243.61
3.714
E
-
17
15
−270.82
5.658
E
-
29Summary of population genetic analysis of RAF nuclear genotyping. (a) Venn diagram illustrating gene flow occurring between populations. The 11Junonia populations are distributed into 3 clusters which show extensive genetic exchange within a cluster, but little or no genetic exchange between clusters. (b) CONTML analysis of the RAF allele frequency data for populations identified in Structure. Clusters identified in (a) tend to be adjacent to one another on the tree, which shows a very strong influence on larval host plant use. Most of the individuals in RAF Populations 2A, 2B, 6A, 6B, and 6C are from species that utilize larval host plants in the order Scrophulariales. Virtually all of theJunonia in the remaining RAF populations are from species that use larval host plants in the order Lamiales.
(a)
(b)Overall, there are 11 populations established by Structure analysis of RAF alleles (Figure5(a)). RAF Population 1 (yellow) primarily includes specimens ofJ. litoralis, but 2J. zonalis, 1J. genoveva, and 1J. evarete specimens also show genetic influence from this population.Population 2A (dark blue) includes all specimens ofJ. wahlbergi andJ. divaricata and all but one of theJ. evarete specimens. Population 2A also includes individuals from all of the other French Guianan and French AntilleanJunonia species. Most of the specimens in Population 2B (light blue) also show genetic influence from Population 2A, but there are 2 CaribbeanJ. zonalis specimens whose primary genetic influence is population 2B. Population 3A (dark purple) is comprised ofJ. genovevaandJ. zonalis specimens. Populations 3B (light purple) and 3C (pink) both includeJ. genoveva,J. neildi, andJ. zonalisspecimens. Population 4 (red) contains onlyJ. genoveva specimens. Population 5 (orange) is primarily composed ofJ. genovevaspecimens, although there is 1J. zonalis specimen that also shows influence from this population. Populations 6A-C are the same as the 3 ArgentinianJunonia 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 (Figure5(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 ofJ. zonalis andJ. neildi from the Caribbean andJ. genovevaandJ. litoralis from French Guiana comprise the remainder of this cluster (along with 1 specimen ofJ. 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 fromJ. evarete,J. wahlbergi, andJ. divaricata from French Guiana, although it also includes specimens from the other 2 French Guianan species, specimens ofJ. neildiandJ. zonalis from the Caribbean, andJ.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 (Figure5(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 withCOI andwingless 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 ofJunonia 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 ofJunonia 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 inJunonia [23, 24, 54]. CertainCOI 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 ofJunonia have been able to provide corroboration supporting the existence of cryptic species associated withCOI 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 ofJunonia is apparently absent from a region [55],COI barcoding should probably be abandoned as a method for distinguishing taxa among New WorldJunonia. At the same time, the segregation of haplotypes A and B, separated by approximately 4% sequence divergence, inJunonia 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 ofCOI 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 CaribbeanJunonia 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 putativeJunonia ring species since several phenotypically distinct forms ofJunonia coexist there [34, 92]. Our results suggest that if there is a region of secondary contact between two ends of a CaribbeanJunonia 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 inJunonia populations.Alleles from the nuclear locuswingless, which were helpful in distinguishing betweenJunonia from Argentina [24], were less effective in this study. This is probably due, at least in part, to the small 137 bp sequence fragment ofwingless that we were able to recover from Caribbean and French GuiananJunonia (versus 402 bp previously recovered for ArgentinianJunonia [24]). While thewingless fragment that was recovered is the most variable portion of the New WorldJunonia 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 ofJunonia from French Guiana and the Caribbean with better-preserved nuclear DNA so that larger portions of thewingless gene could be analyzed. Another factor that may makewingless sequences less useful in French Guiana and the Caribbean is the large number of forms ofJunonia coexisting and possibly interbreeding in this region [23]. Ifwingless coding sequences or sequences closely linked towingless 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 ofJunonia [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 [10–12, 14] may permit us to identify the specific mutations responsible for phenotypic evolution inJunonia 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 ofJunonia 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, ArgentinianJunonia have been conventionally referred to asJ. genoveva hilaris (C. & R. Felder), the light buckeye butterfly, andJ. 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 theJunonia 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 toJ. genoveva and either a genetically disparate population of the same species or an undescribed cryptic species while the dark-coloured population corresponds toJ. wahlbergi. The results of this study suggest that there is no close affinity between the Argentinian forms andJ. genoveva orJ. wahlbergi from French Guiana, so the tentatively assigned scientific names should be revised. We refer to the 2 light-coloured forms asJ. “hilaris” andJ.sp. affin “hilaris” and the dark-coloured form asJ. “flirtea.”The samples from French Guiana and the Caribbean are divided into two major clusters (Figure5). 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 ofJunonia whose larval host plants are in order Lamiales:J. litoralis from French Guiana andJ. neildifrom the Caribbean (both use larval host black mangrove,Avicennia germinans),J. genoveva from French Guiana (larval hostHyptis atrorubens), andJ. zonalisfrom the Caribbean (larval hostsStachytarpheta jamaicensis,S. urticifolia, andLippia nodiflora) [23, 99]. Since allJunonia 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, andJ. wahlbergi.Junonia divaricata andJ. evarete useUtricularia hispida as their larval host, whileJ. wahlbergi andJ. evarete useAgalinis hispidula [23].Junonia coenia from North America feeds on a wide variety of larval hosts in the order Scrophulariales including severalAgalinis 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 WorldJunonia. However, it should also be noted that there are exceptions to this congruence: 1 specimen ofJ. evarete clustered with the Lamiales feeders while 12 specimens ofJ. genoveva, 5 specimens ofJ. zonalis, and 5 specimens ofJ. neildi clustered with the Scrophulariales feeders (Figure 5). This might be interpreted as evidence for some gene flow between these two lineages.Some species ofJunonia were not readily distinguishable from one another in the Structure analysis of the RAF data.Junonia coenia,J. divaricata,J. evarete, andJ. 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 theseJunonia species will be possible. While we were not able to reliably separateJ. zonalis andJ. 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, 26–28].The other pattern emerging from the analysis of RAF genotypes is for individuals of one species to be spread across multiple RAF populations (Figure5).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, otherJ. 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 whileJ. 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 ofJunonia existing simultaneously in Buenos Aires, Argentina, may be due to mass migrations of individuals from geographically disparate areas [24]. Mass migration ofJunonia is a phenomenon that has been documented in Argentina and elsewhere in the New World [100, 101]. However, mass migration ofJunonia 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 namedJunonia taxa from French Guiana and the French Antilles have been described relatively recently [23, 47–49], 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 WorldJunonia tested in lab crosses are interfertile and produce viable fertile hybrids [39–41], 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 sympatricJunonia 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 inJunonia. Several species ofJunonia have been tested forWolbachia (includingJ. evarete from Panama), but thus far infections have only been detected in AsianJ. 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 inJunonia at that locality unless there is assortative mating between forms.In addition to explaining patterns of genetic overlap amongJunonia RAF populations, we have to explain the possible subdivision of someJunonia 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 ofJunonia 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 [109–113]. In some cases, it has been suggested that this has been a driver for reproductive isolation and incipient speciation in many insects [113–117].Most New WorldJunonia 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 ofJunonia 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, femaleJ. 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 ofJunonia exist, this may explain some of the extensive population structure seen in some species such asJ. genoveva. New WorldJunonia host plants contain iridoid glycoside secondary compounds, the presence of which may be a necessary precondition for use as a host byJunonia [120]. This may help identify possible additional larval host plants for South American forms ofJunonia.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 multipleJunonia taxa occur, one of the authors (Jeffrey M. Marcus) has observed differences in habitat use. In mangrove swamps along the coast of Florida, USA, mostJ. coenia males appear to patrol mating territories in clearings without trees or other vertical habitat structure. In contrast, males of a form that resemblesJ. 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 pigmentedJ. “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 AmericanJunonia 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 differentJunonia 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 cooccurringJunonia 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 ofJunonia, but there are seasonal differences that may contribute to an allochronic mechanism for persistence. In French Guiana, the foliage ofJunonia 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 ofJ. divaricata,J. evarete, andJ. wahlbergi coincide temporally with each other and with the presence of larval hosts in the Scrophulariales while the flight times ofJ. genoveva,J. litoralis,J. neildi, andJ. zonalis are less seasonally restricted [23]. This difference in phenology may contribute to the continued distinctiveness of the two majorJunonia 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 aboutJunonia 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 ofJunonia [26, 108]. Variation in courtship flight patterns in otherJunonia forms has not yet been documented. There are also differences in the colour patterns of New WorldJunonia 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 ofJunonia 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. WhileCOIbarcodes are of limited utility, nuclearwingless sequences and RAF genotyping are effective at identifying some individual species ofJunonia and have been very helpful in examining the relationships ofJunonia 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 twoJunonia 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 WorldJunonia 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.



---

*Source: 897596-2014-10-12.xml*