Abstract
Cytochrome c oxidase subunit I (COI) sequences were used to estimate demographic histories of populations of the buckeye butterfly Junonia genoveva (Cramer) from Costa Rica and Mexico. Previous studies have revealed significant structure between populations of J. genoveva from coastal regions of northwestern Mexico, which utilize black mangrove Avicennia germinans (Acanthaceae) as a larval host plant, and inland populations from Costa Rica that feed on different hosts in the families Acanthaceae and Verbenaceae. The Mexico population of J. genoveva reported on here is located near the Northern limit of black mangrove habitat on the Pacific coast of North America and is hypothesized to have been established by northward migrations and colonization from southern source populations. The mismatch distribution, Bayesian skyline analyses, and maximum likelihood analyses carried out in FLUCTUATE were used to estimate changes in female effective population size () over time in the two populations. Differences found in COI haplotype diversity, present-day , and the timing of population expansions are consistent with the hypothesis that the Mexico population of J. genoveva is the more recently evolved.
1. Introduction
The genus Junonia Hübner comprises approximately 30+ species of butterflies commonly known as buckeyes and pansies which are distributed worldwide, predominately in the tropics. The New World fauna is thought to have originated about 2–4 million years ago (Ma) from colonizing individuals from Africa or Asia [1]. The number of species of Junonia in the New World is uncertain, but at least three, including J. evarete (Cramer), J. genoveva (Cramer), and J. coenia Hübner, are reported for North America [2]. Junonia evarete and J. genoveva are also found in South America and the Caribbean.
Molecular studies using both mitochondrial and nuclear DNA suggest that two distinct genetic lineages (clades) of Junonia began to differentiate in the New World shortly after colonization [1, 3]. Evidence for the relatively recent speciation, together with the pronounced intraspecific phenotypic variability and tendency to hybridize [4], has caused much taxonomic confusion within the genus, especially in distinguishing J. evarete and J. genoveva [5]. The populations of Junonia treated here typically have been assigned to J. evarete [6–8], but because these populations are genetically distinct from J. evarete, it has been suggested that they be provisionally reassigned to J. genoveva [3]. Significant structure was found between a coastal population of J. genoveva from northwestern Mexico and an inland population from Costa Rica () [3], possibly owing to the large geographic separation (approximately 3250 km; Figure 1(a)) and different host plant preferences. The northwest Mexico population utilizes black mangrove, Avicennia germinans (L.) L. (Acanthaceae), as a larval host [9], whereas the Costa Rica population feeds on Dyschoriste valeriana Leonard (Acanthaceae) and Stachytarpheta jamaicensis (L.) Vahl (Verbenaceae) [10]. The higher COI haplotype and nucleotide diversities found in the Costa Rica population compared with the Mexico population suggested a scenario in which dispersal and colonization of J. genoveva proceeded northward into Mexico from populations originating in Central or South America [3]. Here, we use several tests of COI sequence data from the two populations to obtain estimates of demographic histories in an attempt to provide further insight into the northward expansion hypothesis and the present-day distribution of J. genoveva in Mexico.
(a)
(b)
2. Materials and Methods
2.1. Sampling
Adults of J. genoveva from northwestern Mexico were collected at Estero del Soldado, near Guaymas, Sonora, Mexico (Figure 1(b)), and at nearby San Carlos [9]. One adult was reared on black mangrove from a first instar larva obtained at Estero del Soldado. Sample size was . Adults from the Costa Rica population (), listed as J. evarete [10], were obtained from larvae collected from the Area de Conservación Guanacaste (ACG), Guanacaste province in northwestern Costa Rica (Figure 1(a)) and reared on Dyschoriste valeriana or Stachytarpheta jamaicensis.
To estimate seasonal abundance of J. genoveva at Estero del Soldado, we conducted a weekly survey of adults from February 2011 to February 2012. Adults encountered on a 0.5 km dirt road adjacent to the mangrove forest (Figure 1(b)) were averaged over two-week intervals. The pattern of abundance (Figure 2) is in general agreement with Brown et al. [7] who found that J. genoveva (as J. evarete) in Baja California Sur, Mexico flies from September through February, but is most common from September through November. Although we cannot state with certainty the number of generations per year in J. genoveva, our data are consistent with multiple broods. For the demographic tests, we have chosen a conservative value of two generations per year for the Guaymas population and have assumed the same value for the Costa Rica population (see Section 4).
2.2. Molecular Protocol
Details on standard procedures used for extraction of total genomic DNA from butterfly legs and amplification of the COI gene segment are found elsewhere [3, 12]. The 658 bp segment corresponds to the COI barcode region [13]. GenBank accession numbers are JQ430692–JQ430710 (Guaymas, Mexico) and GU334034, GU334037, GU157280–GU157292, and GU157294–GU157300 (ACG, Costa Rica). Specimen voucher codes and geographic coordinates for corresponding GenBank accession numbers are found in Pfeiler et al. [3]. Detailed information on Costa Rica specimens can be obtained by entering voucher codes into the ACG database [10].
Twelve additional COI barcode sequences were available in GenBank for Junonia from four states in southern Mexico: Morelos (JQ430731, JQ430733), Campeche (GU659427, GU659432, GU659436), Yucatán (GU659429–GU659431), and Quintana Roo (HQ990188, GU659425, GU659426, GU659435). Although these sequences are currently assigned to J. evarete, our analyses showed that they belonged to the J. genoveva genetic lineage. These samples were not used in the demographic analyses owing to low sample sizes from each region, but their geographic distribution (Figure 1(a)) provides additional support for our colonization hypothesis.
2.3. Data Analysis
Calculations of genetic diversity indices and Tajima’s [14] D were performed in DnaSP version 5.00.04 [15]. Fu's [16] Fs neutrality tests were conducted in ARLEQUIN version 3.5.1.3 [17] using 1000 simulations. In addition to assessing whether nucleotide polymorphisms deviate from expectation under neutral theory, Fu's Fs test is also useful for detecting signatures of population expansions, which lead to large negative values in the test statistic [16, 18]. The significance of Fs at the 0.05 level was indicated when P values were <0.02 [17].
Estimates of changes in effective female population size () over time for J. genoveva were obtained in FLUCTUATE version 1.4 [19] and from Bayesian skyline analysis implemented in BEAST version 1.2 [20]. FLUCTUATE provides simultaneous maximum likelihood estimates of the mutation parameter (θ) and the exponential population growth parameter (g), where , and μ is the neutral mutation rate per site per generation. To estimate μ, we assumed 2.3% pairwise sequence divergence per million years for the COI gene [21] and two generations per year, yielding a single lineage value of . Details on the program settings used in FLUCTUATE are given in Pfeiler et al. [22]. Bayesian skyline analysis utilizes Markov chain Monte Carlo (MCMC) sampling of sequence data to estimate a posterior distribution of effective population size through time [20]. Bayesian skyline analyses were run under the conditions of the HKY + Γ model (four gamma categories) using the mean mutation rate per site per generation (μ) described above. Ten million iterations of the MCMC chains were run, sampling every 1000 iterations. The Bayesian skyline plots were generated with TRACER version 1.2.1 [20]. To provide confirmation of population expansions detected in FLUCTUATE and Bayesian skyline analyses, demographic histories of both populations were also tested by analyzing the distribution of pairwise sequence differences, also known as the mismatch distribution [23, 24], using ARLEQUIN. The significance of the estimated parameters of the sudden expansion model of the mismatch distribution is obtained from the sum of square deviations (SSD) statistic and the raggedness statistic (rg), and their corresponding P values. The sudden expansion model is rejected when P is <0.05.
3. Results
Genetic diversity indices for the 658 bp COI segment reported previously [3] revealed that haplotype diversity (h) and nucleotide diversity (π) were higher in the Costa Rica population (; ; ) than in the Guaymas, Mexico population (; ; ). Tajima's D values also were not significant in either populations. In the present study, a significant negative value for Fu's Fs was found for the Costa Rica population (; ). Although the value for the Guaymas population was not significant (; ), it approached the cutoff level for significance. The results suggest a historical population expansion for the Costa Rica population, which as we show below is consistent with results obtained from the other demographic tests. These tests also suggest a historical population expansion for the Guaymas population, but in all cases the evidence for this expansion is weaker than for the Costa Rica population.
The results from FLUCTUATE (Table 1) showed a positive exponential growth parameter (g) that was significantly different from zero in both populations of J. genoveva, indicating an increase in population size, although the standard deviation of g was much higher in the Guaymas population. Also, effective female population size () was much larger in J. genoveva from Costa Rica than in the population from Guaymas (Table 1).
Bayesian skyline analyses (Figure 3) also indicated that both populations of J. genoveva have increased in size, with the present day again much larger in the Costa Rica population (Figure 3(a)). The timing of the expansion of the Guaymas population (Figure 3(b)) was estimated to be more recent (~86,000 years before present, BP) than the Costa Rica population (~400,000 BP).
(a)
(b)
Results from the mismatch distribution are shown in Table 2 and Figure 4. For the Costa Rica population (Figure 4(a)), the plot of the observed distribution of pairwise differences among COI haplotypes was unimodal and agreed well with the expected distribution for a population that has undergone an expansion. A unimodal curve also was seen in the Guaymas population, but the fit with the expected distribution was not as good (Figure 4(b)). For both populations, the test statistics SSD and rg were small and not statistically significant (Table 2), indicating that the sudden expansion model could not be rejected for either population. The P values for SSD and rg for the Guaymas population, however, were near the 0.05 cutoff level for significance (Table 2), suggesting that evidence for population size increase is weaker in this population, consistent with the results from FLUCTUATE.
(a)
(b)
Based on the results from the mismatch distribution (Table 2), we estimated the timing of the population expansion in both the Costa Rica and Guaymas populations of J. genoveva using the equation , where is a moment estimator of mutational time, t is the number of generations since the expansion, and u is the mutation rate for the entire gene segment of 658 bp [23]. Assuming 2.3% pairwise divergence per million years [21], the mean mutation rate per site per generation in the 658 bp segment for a single lineage is equal to or . The estimated time to the population expansion in the Costa Rica population (with 95% confidence intervals) was 132,450 (52,318–375,500) generations ago, and was 70,199 (13,245–145,030) generations ago in the Guaymas population. Although the confidence intervals are large, these calculations reveal the same trend as seen in Bayesian skyline analysis in which the timing of population expansion in the Costa Rica population predates that of the Guaymas population.
4. Discussion
Our results show that female effective population size () and genetic diversity indices are larger in the Costa Rica population of J. genoveva compared to the northwestern Mexico population at Guaymas, consistent with the hypothesis that the Guaymas population was founded by migrating individuals from Central America as depicted in Figure 1(a). In addition, the presence of J. genoveva throughout the Yucatán Peninsula and the state of Morelos revealed by barcode analysis of available GenBank sequences (listed as J. evarete) suggests that colonization from Central America source populations was widespread in Mexico. Additional sampling will be required to determine if the J. genoveva lineage also colonized northeastern Mexico and southern USA, but this scenario seems probable based on the data presented here.
Population expansions in the Guaymas and Costa Rica populations of J. genoveva both dated to the Pleistocene, but differences in the timing of the expansions (Figure 3) are also consistent with the conclusion that the Guaymas population is the more recent. These differences, however, are only rough estimates that depend on the assumption of a standard 2.3% molecular clock and two generations per year in both populations. Different numbers of generations per year would result in changes in the estimated timing of the expansions, but it is unlikely that our conclusions would be affected. For example, if we assume a value of three generations per year in the Costa Rica population, the estimated time of the expansion would increase by about 50% to ~600,000 years BP, still within the Pleistocene and much earlier than the population expansion at Guaymas (~86,000 BP). In an unlikely scenario of one generation per year in the Costa Rica population and three generations per year at Guaymas, the Costa Rica population expansion would still predate the Guaymas expansion (~200,000 and ~130,000 years BP, resp.).
The apparent absence of J. evarete in Mexico and Central America as suggested by barcode analysis is noteworthy, as this species is reported to be a common inhabitant of the region [6–8]. More extensive sampling may ultimately reveal J. evarete, but all COI barcodes obtained to date (), from samples covering a broad geographic area from northwestern Mexico to Panama, and including both coastal and inland populations utilizing different larval host plants, clearly belong to the J. genoveva lineage ([3], present study).
Specimens of Junonia from Baja California Sur listed as J. evarete [7, 8] show an ecological association with the mangrove habitat, in addition to showing similar patterns in wing color and seasonal abundance compared with the Guaymas population, leading to the suggestion that they should be reassigned to J. genoveva [3, 9]. Given this evidence, we show colonization of Baja California Sur by J. genoveva on Figure 1(a), although no COI barcodes are available for confirmation. Our demographic evidence suggesting a late Pleistocene origin for the mainland population at Guaymas would also imply that putative colonization of J. genoveva in Baja California Sur occurred by over-water dispersal from the mainland, as the separation of the Cape Region of the Baja California peninsula from the mainland is thought to have occurred during the Pliocene, about 3-4 Ma [25]. In support of this hypothesis, J. villida (Fabricius), a species closely related to the New World Junonia [1], is known to be a strong over-water disperser, showing a broad geographic distribution among Pacific islands [26, 27].
Acknowledgments
The authors thank T. Hernández Mendoza and M. Polihronakis-Richmond for help with this project. They are especially grateful to D. H. Janzen, W. Hallwachs and N. Wahlberg for kindly sharing their unpublished data. They also thank R. A. Bailowitz, C. Brévignon, J. Calhoun, and K. Hansen for helpful discussions and providing insights into the taxonomy of Junonia. This research was supported by NSF Grant DEB-0346773 to T. A. Markow.