Independent divergence of 13- and 17-y life cycles among three periodical cicada lineages
Edited by Mary Jane West-Eberhard, Smithsonian Tropical Research Institute, Ciudad Universitaria, Costa Rica, and approved February 22, 2013 (received for review November 16, 2012)
Commentary
April 9, 2013
Abstract
The evolution of 13- and 17-y periodical cicadas (Magicicada) is enigmatic because at any given location, up to three distinct species groups (Decim, Cassini, Decula) with synchronized life cycles are involved. Each species group is divided into one 13- and one 17-y species with the exception of the Decim group, which contains two 13-y species—13-y species are Magicicada tredecim, Magicicada neotredecim, Magicicada tredecassini, and Magicicada tredecula; and 17-y species are Magicicada septendecim, Magicicada cassini, and Magicicada septendecula. Here we show that the divergence leading to the present 13- and 17-y populations differs considerably among the species groups despite the fact that each group exhibits strikingly similar phylogeographic patterning. The earliest divergence of extant lineages occurred ∼4 Mya with one branch forming the Decim species group and the other subsequently splitting 2.5 Mya to form the Cassini and Decula species groups. The earliest split of extant lineages into 13- and 17-y life cycles occurred in the Decim lineage 0.5 Mya. All three species groups experienced at least one episode of life cycle divergence since the last glacial maximum. We hypothesize that despite independent origins, the three species groups achieved their current overlapping distributions because life-cycle synchronization of invading congeners to a dominant resident population enabled escape from predation and population persistence. The repeated life-cycle divergences supported by our data suggest the presence of a common genetic basis for the two life cycles in the three species groups.
Sign up for PNAS alerts.
Get alerts for new articles, or get an alert when an article is cited.
Periodical cicadas (Magicicada) in the eastern United States represent one of the most spectacular life history and population phenomena in nature (1–10). These periodical cicadas spend most of their lives (13 y in the south, 17 y in the north) as underground juveniles except for a brief 2- to 4-wk period when adults emerge simultaneously in massive numbers. With few exceptions, at any given location, all of the periodical cicadas share the same life cycle and emerge on the same schedule, forming a single-year class referred to as a “brood.” Surprisingly, each brood consists of multiple species from three species groups (Decim, Cassini, Decula).
These three groups were considered to have diverged from each other allopatrically and to have later become sympatric and formed 13- and 17-y life cycles (2). The prolonged, prime-numbered life cycles were hypothesized to have evolved in response to Pleistocene climatic cooling (9, 11) to avoid the adverse effect of low population density on mating success (9, 12, 13). Another view hypothesized that the long synchronized life cycles evolved in association with the predator avoidance strategy (2, 4, 8) and that this took place before both the glacial periods and the split of the three species groups (10) based on approximate genetic distances among species groups (8). To test these hypotheses, phylogenetic information about the relationships of species, broods, and populations is essential. However, phylogenetic studies of Magicicada have been largely limited to the Decim group (6, 14–18), and only rough divergence times among species have been inferred (8, 10). Thus, a comprehensive molecular phylogeny covering all of the extant broods, their phylogeography, and divergence time has been lacking until now.
We conducted molecular phylogenetic and population genetic analyses by using nuclear and mitochondrial DNA markers for samples collected over a 30-y period (1978–2008). These samples represent all 15 extant broods and all known species (Table S1).
Results and Discussion
Phylogenetic analyses and divergence time estimation based on four nuclear and three mitochondrial genes (Table S2) clearly reveal the monophyly of each of the three species groups (Decim, Cassini, and Decula) and the sister relationship between Cassini and Decula (Fig. 1). A Bayesian relaxed clock analysis shows that the three species groups diverged 3.9 Mya. Initially the Decim group diverged from the ancestor of Cassini + Decula, then Cassini and Decula separated 2.5 Mya (Fig. 1). The mitochondrial gene genealogy further shows divergence associated with regions, and partly with life cycles (Fig. 2, Figs. S1 and S2, and Table S3). We distinguished four mitochondrial haplotype groups in Decim and three each in Cassini and Decula (Fig. 2 A–C). The geographic distributions of mitochondrial haplotype groups within each species group show similar divisions among eastern, middle, and western regions (Fig. 2 D–F). In Decim, there is also a major divergence between northern and southern groups corresponding to formally distinguished mitochondrial lineages A and B, respectively (15). Group A is divided into three groups, Ae, Am, and Aw, which occur in the eastern, middle, and western parts of the United States east of the Great Plains, respectively. Populations of Ae and Am exhibit a 17-y cycle (Magicicada septendecim); those of Aw both 17-y (M. septendecim) and 13-y (Magicicada neotredecim) cycles; and those of B show only a 13-y cycle (Magicicada tredecim). At the boundary of Ae and Am, a few populations possess both Ae and Am haplotypes. The Cassini group consists of three haplotype groups, Ce, Cm, and Cw, again occurring in the eastern, middle, and western regions, respectively (Fig. 2). Populations with Ce and Cm show 17-y cycles only (Magicicada cassini), whereas those of Cw show both 17- and 13-y cycles (M. cassini and Magicicada tredecassini, respectively). Lastly, the Decula group shows the least divergence of mitochondrial haplotypes, but yet there are differences in haplotype among the eastern (De), middle (Dm), and western (Dw) regions (Fig. 2). Each haplotype group is associated with both 17- and 13-y cycles (Magicicada septendecula and Magicicada tredecula, respectively), suggesting that life cycle divergence occurred independently in the three regions.
Fig. 1.
Fig. 2.
The boundaries of haplotype groups roughly coincide among the three species groups (Fig. 2 D–F). In each species group, an isolation-by-distance pattern was detected by a Mantel test (Decim, Mantel’s r = 0.223, P = 0.0007; Cassini, r = 0.349, P = 0.0001; Decula, r = 0.377, P = 0.0001), suggesting that contiguous range expansion resulted in the present distribution. Parts of broods of each species group were associated with two or three haplotype groups and occurred in different geographic regions, suggesting multiple origins of the same broods even within a single species (Fig. 2 and Fig. S2). In Decim, broods VI, IX, X, and XIV are found in Ae and Am, and broods XIX and XXIII are found in Aw and B. In Cassini, broods V and IX are found in Ce and Cm, and broods X and XIV are in Ce, Cm, and Cw. In Decula, broods V and XIV are found in De and Dm, broods VI and XXIII in De and Dw, and broods XIX in De, Dm, and Dw. Thus, the broods of multiple lineages are shared among the three species groups and show similar geographic patterns. The divergence pattern of the two life cycles differs among the three species groups in relation to haplotype groups (Fig. 2). Although genetic differentiation between 13-y M. tredecim and 17-y M. septendecim in Decim are evident, naturally, because these belong to different haplotype groups [Table S4; analyses of molecular variance (AMOVA)], significant differentiation between 13-y M. neotredecim and 17-y M. septendecim within the haplotype group Aw was also detected by an AMOVA (Table S4; FCT = 0.234; df = 1, P < 0.001). However, there was not significant differentiation between 13- and 17-y cicadas in Cassini haplotype group Cw or in the whole Decula group (Table S4; Cassini Cw: FCT = 0.068, df = 1, P > 0.05; Decula: FCT = −0.039, df = 1, P > 0.05).
The analysis of amplified fragment length polymorphism (AFLP) markers show results generally consistent with the mitochondrial data results (Fig. 3). Each species group showed a significant differentiation in terms of the fixation index among groups defined by life cycle and regions associated with mitochondrial haplotype groups (Decim, Fst = 0.0241, P = 0.0000; Cassini, Fst = 0.0099, P = 0.0000; Decula, Fst = 0.0109, P = 0.0040). In Decim, the 13-y M. tredecim (haplotype group B) was significantly differentiated from other groups as indicated by pairwise Fst values, whereas another 13-y group M. neotredecim (haplotype group Aw) and 17-y M. septendecim groups were not differentiated from one another except for one case (Fig. 3). In Cassini, the closest relationship between 13-y M. tredecassini and 17-y M. cassini within the haplotype group Cw was revealed (Fig. 3). In Decula, the six groups divided by life cycle and haplotype group showed relationships different from those inferred by using mitochondrial data; 13-y M. tredecula in haplotype groups De and Dm were closely related to each other, as were 17-y M. septendecula in these haplotype groups (Fig. 3). Note, however, that haplotype groups De and Dm show only 1-bp difference in the mitochondrial sequence (Fig. 2C).
Fig. 3.
The estimated age using a Bayesian relaxed clock analysis of the most recent common ancestor (MRCA) of all Decim mitochondrial haplotypes is 0.53 Mya, after which time it separated into clades A (Ae+Am+Aw) and B (Fig. 1). The MRCA of each haplotype group is more recent (0.16–0.08 Mya). The MRCA of the Decim Aw group containing 13-y (M. neotredecim) and 17-y (M. septendecim) cicadas is 0.12 Mya (Fig. 1). In the Cassini group, the MRCA of all mitochondrial haplotype groups is traced back to 0.32 Mya, and the MRCA of Cw back to 0.16 Mya. Lastly, mitochondrial haplotypes of the Decula group had their MRCA 0.23 Mya. The divergence time of 13- and 17-y cicadas within the same mitochondrial lineage was estimated together with the gene flow between them using the isolation-with-migration approach (Fig. 4 and Fig. S3). In Cassini Cw, divergence time was ∼23 kya, whereas in Decim Aw, divergence time was ∼10 kya. In Decula, the divergence time was not clearly estimated (close to zero) due to the lack of haplotype differences between M. septendecula and M. tredecula, so their life cycle divergence must have occurred quite recently.
Fig. 4.
The Bayesian skyline plot (BSP) based on mitochondrial gene sequence data (Fig. 5) as well as that based on both mitochondrial and nuclear gene sequence data (Fig. S4) revealed that the population sizes of both the Decim and the Cassini groups were relatively small during the last glacial period, probably due to a population bottleneck, but increased markedly after 10,000 y ago following the abrupt rise of temperature in the Holocene (19). By contrast, there was no increase in population size for the Decula group. This difference in the demographic history corresponds with the fact that except in some southern populations, today Decula is generally rare compared with Decim and Cassini (8).
Fig. 5.
Our results are broadly consistent with the previous idea that an ancestor of all Magicicada diverged into three species allopatrically, and later, the three became sympatric and each species independently diverged into 13- and 17-y cicadas (2). Surprisingly, however, the divergence of 13- and 17-y cicadas was asynchronous among the species groups and occurred repeatedly even within a species group. This finding is all of the more interesting given that each species group shows similar eastern, middle, and western phylogeographic divisions (Fig. 2) similar to post-Pleistocene patterns observed in other North American taxa (20), suggesting that the three Magicicacda groups shared multiple refugia during the last glacial maximum. The Decim group showed the deepest (0.53 Mya; middle Pleistocene) divergence in extant haplotype groups between north and south (i.e., Ae+Am+Aw vs. B). Meanwhile, the divergence between 13- and 17-y cicadas in Decim Aw, Cassini Cw, and Decula occurred during or after the last glacial maximum (<23 kya), possibly following population bottlenecks, suggesting that the life cycle divergence in Magicicada is closely associated with global climatic fluctuations and shorter growing seasons in the north versus the south. The divergences into the present 13- and 17-y cicadas have occurred since the middle Pleistocene when climatic fluctuation became prominent (19), consistent with the glacial-period-origin hypothesis of the life cycles (9, 11). However, the splits into species groups date to the Pliocene, suggesting that the origin of the life cyles was not directly related to the onset of the Pleistocene glacial period given our reasoning below (10).
We suggest that the recent, repeated life cycle divergences are caused by shifts between the alternative life cycles. We further suggest that the genetic basis of these life cycles and shifts is the same in each species group because it is unlikely that the same regulatory system evolved three times. The origin of M. neotredecim has been proposed to be a recent life cycle shift of M. septendecim from a 17-y to a 13-y cycle by genetic assimilation of life cycle plasticity, perhaps facilitated by other 13-y cicada “nurse broods” that protected them from predation (4, 17, 21). The plasticity of Magicicada life cycle length has been suggested based on records of off-schedule emergences (22). Supporting the nurse brood model is the fact that M. neotredecim is found in two of three extant 13-y cicada broods, XIX and XXIII with broad ranges, but not in XXII with a narrow range far away from the range of M. neotredecim. Thus, the origin of M. neotredecim may have involved some influence by preexisting 13-y broods XIX and XXIII. Introgressive hybridization between M. septendecim and M. tredecim (23, 24) has also been proposed to explain the origin of M. neotredecim. Our results corroborate the previously reported genetic closeness of M. neotredecim to M. septendecim and not to M. tredecim (15–18), more consistent with the plasticity/nurse-brood model (but see ref. 25 for the theoretical possibility of the introgressive hybridization model). In the Cassini and the Decula groups, the hybridization hypothesis cannot be applied because there is no evidence of secondary contact between cicadas with different life cycles of the same species group. These species group may have been monotypic (either 13- or 17-y) at the beginning of the last glacial period and may have changed their life cycles to match those of preexisting nurse-brood populations (e.g., of Decim) with which they became sympatric. Natural selection would have promoted synchronization of invading populations to resident populations because invaders would gain protection from predation (2, 4, 8) and, consequently, avoid Allee effects (failure to reproduce due to low population density) (9). The repeated shifts between the two prime-numbered life cycles in all three Magicicada species groups suggests a common genetic basis that evolved before the origin of the species groups, 3.9 Mya. Our results support the idea that life cycle plasticity has been a creative force in the evolution of Magicicada (26). However, the regulatory system of Magicicada life cycles and the mechanism of putative life cycle shift are completely unknown and require future comparative genomic studies of 13- and 17-y periodical cicadas.
Materials and Methods
Sampling.
Adults of Magicicada were collected from 27 states in the United States from 1978 to 2008 (Table S1) and stored in freezers. Total genomic DNA was extracted from leg muscles of each ethanol-fixed specimen by using the Wizard Genomic DNA Purification Kit (Promega). For Decim, Cassini, and Decula, 332, 238, and 165 individuals, respectively, were used. Specimens are kept in the Department of Zoology, Kyoto University, as vouchers and for future analyses of DNA. We used the sister genera of Magicicada (27) as outgroup taxa, specifically, Tryella crassa, Tryella graminea, Tryella burnsi, Aleeta curvicosta, and Aleeta curvicosta.
Gene Sequence Analysis.
We sequenced partial gene regions of nuclear 18S rRNA (18S; 808-bp sequence), wingless (Wg; 376-bp exon sequence), and Elongation Factor 1-alpha (EF1-a exon and intron sequence; 1,302–1,303-bp for Magicicada; 1,270–1,616-bp for the outgroup) and Calmodulin (Cal; 385–394-bp intron sequence) gene and a 1,832-bp mitochondrial DNA region encompassing partial sequences of cytochrome oxidase subunit I and II genes (COI, COII) and a sequence of tRNA (Leu) in between. Primers used for PCR amplification and direct sequencing are given in Table S2. Direct sequencing of PCR products used an ABI3130xl sequencer (Applied Biosystems). When sequence data for the COI-II sequence were affected by putative nuclear mitochondrial DNA (numt), we used specific internal primers (Table S2). Samples of each species group representing different broods and regions were sequenced for all gene regions. The samples consisted of 25 from the Decim group, 18 from the Cassini group, and 19 from the Decula group. Two outgroup specimens (T. crassa and A. curvicosta) were also sequenced. For the mitochondrial DNA, sequences were obtained for 329, 238, and 165 specimens of Decim, Cassini, and Decula, respectively (Table S1). All nuclear gene sequences and all mitochondrial haplotype sequences detected have been deposited in DNA Data Bank of Japan database (accession nos.: 18S, AB740543–AB740606; Cal, AB740607–AB740670; EF1-a, AB740671–AB740734; Wg, AB740735–AB740798; COI-II, AB740799–AB740917; see also Table S3 for COI-II haplotypes).
Phylogenetic Analysis.
We conducted a simultaneous analysis of five gene sequences with 64 operational taxonomic units (OTUs) using an optimal data-partitioning scheme (28). The nuclear gene sequence data were divided into nine subsets: 18S and Cal (intron); four of EF1-a (first, second, and third position, and intron); and three codon positions of Wg. Mitochondrial gene sequence data were divided into seven subsets: three codon positions of COI; tRNA; and three codon positions of COII. These 16 subsets were used in a heuristic search to identify the best partitioning scheme and evolutionary model for each partition according to the Bayesian information criterion by using the program PartitionFinder version 1.0.1 (28). The best scheme consisted of six partitions: (i) 18S + EF1-a second + Wg second [Substitution model: JC (Jukes-Cantor)]; (ii) COI first + tRNA(Leu) + COII first + Cal + EF1-a intron [HKY (Hasegawa-Kishino-Yano) + G]; (iii) Wg third [TVM (transversion model) + G]; (iv) EF1-a first + Wg first [F81 (Felsenstein 1981)]; (v) COI second + COII second + EF1-a third (HKY+I); and (vi) COI third + COII third (HKY+I). For mitochondrial gene sequences, a larger data set with 119 OTUs (114 unique sequences of Magicicada and 5 outgroup sequences) was analyzed. The optimal partitioning scheme consisted of three partitions: (i) COI first + tRNA(Leu) + COII first [TrN (Tamura-Nei) + G]; (ii) COI second + COII second (TrN+I); and (iii) COI third + COII third (TrN+I). For each sequence dataset, a partitioned maximum-likelihood analysis was conducted by using Treefinder version October 2008 (29). Bootstrapping analysis with 1,000 replications was performed also by using Treefinder using the ML tree as the starting tree.
Haplotype Network and Population Genetic Analyses of Mitochondrial Gene Sequences.
The relationships among mitochondrial haplotypes were assessed by constructing statistical parsimony networks with a 95% connection limit using TCS version 1.21 (30). Haplotype groups within each species group were assigned based on monophyly or the segregation from other groups with two or more missing haplotypes. Mantel tests for the isolation-by-distance trend in pairwise genetic distance between populations were conducted by using the R package (31). AMOVA for genetic differentiation by haplotype groups, broods, and life cycles were performed by using Arlequin version 3.11 (32).
Divergence Time Estimation.
Estimation of divergence times was conducted with all nuclear and mitochondrial DNA sequence data by using BEAST version 1.7.2 (33). No fossil or geographic evidence for calibrating nodes is available for the Magicicada phylogeny. Therefore, we attempt to determine a plausible range of age for the node of the MRCA of Magicicada by using a recently proposed evolutionary rate for insect mitochondrial genes (ref. 34; Fig. 3), which was estimated by compiling sequence divergence rates at various time spans based on sequence divergence corrected for rate heterogeneity among sites using a gamma distribution and fitting them to a time-dependent rate equation (35): rate of sequence divergence (%) per million year = 17.256 exp[−1.157t] + 2.0968. We can estimate divergence time for a given sequence divergence with this rate by using equation 7 in ref. 35, which relates sequence divergence to divergence time. For the mitochondrial sequences of Magicicada, the optimal substitution model was GTR+G based on the Akaike's Information Criterion according to MrModeltest version 2.3 (36). To determine the node age for the MRCA of Magicicada, the GTR+G corrected sequence divergence between Decim and Cassini + Decula, 0.2351 (SD = 0.0223), was converted to 4.16 My by using the above time-dependent sequence divergence rate. In divergence time estimation with BEAST, the node age prior was set as a normal distribution function with mean = 4.16 and SD = 0.3943 (proportional to original SD). Ref. 34 also provides another time-dependent clock based on uncorrected P sequence divergence from different datasets. Using this clock equation and an uncorrected sequence divergence of 0.0844 ± 0.0034 (SD) for Magicicada MRCA, we obtain 3.57 ± 0.1442 Mya. Because this alternative estimate is included well in the range (3.4–4.9 Mya for 95% HPD interval) of the previous estimate (assuming a normal distribution with mean = 4.16 and SD = 0.3943), our age prior included the age estimation based on a conservative substitution rate without accounting for rate heterogeneity among sites. In the BEAST analysis, data partitioning and substitution models followed the previous phylogenetic analysis. Both substitution model and clock model were unlinked among partitions, but the tree was linked. We tested runs with the strict clock model and the uncorrelated lognormal relaxed-clock model, and compared these by Bayes factor for marginal likelihood using TRACER version 1.5 after the runs. The BEAST Markov chain Monte Carlo run was conducted for 100 million generations sampling every 10,000th generation. The initial 1 million generations were discarded as burn-in when we obtained the maximum clade credibility tree. Based on Bayes factor, we selected uncorrelated lognormal relaxed-clock model.
AFLP Analysis.
AFLP analysis (37) was performed for 237, 95, and 96 specimens of Decim, Cassini, and Decula, respectively, encompassing all broods from each species group (Table S1). We used a plant mapping kit (Applied Biosystems) and six primer combinations for selective amplification (EcoRI+/MseI+: ACT/CTG, AGG/CAC, AGC/CAT, ACA/CAC, AGG/CTA and ACC/CAT). The amplified fragments were electrophoresed on an ABI 3130xl sequencer and binary-coded by using GeneMapper version 4.0 (Applied Biosystems). To ensure high reliability of analyzed AFLP loci, every specimen was genotyped twice, and a total of 440 loci with >90% repeatability was used in the analysis. The differentiation in AFLP loci among groups defined by life cycle and regions associated with mitochondrial haplotype groups within each species group was assessed by using pairwise Fst using AFLP-surv version 1.0 (38). The significance of Fst was assessed by 1,000 permutations. We controlled the false-positive rate by using the Benjamini and Hochberg method (39) to determine statistical significance of multiple pairwise Fst values between groups. We also obtained pairwise Nei’s distances between groups and 1,000 bootstrap distance matrices by using AFLP-surv. We constructed neighbor-joining trees and a bootstrap consensus tree by using PHYLIP version 3.69 (40).
Divergence and Gene Flow Between 13- and 17-y Cicadas.
We applied the isolation-with-migration approach (41) by using the IMa2 program (42) to explore the divergence time between 13- and 17-y cicadas within the same mitochondrial lineage, i.e., those in Aw of Decim (M. septendecim and M. neotredecim), those in Cw of Cassini (M. cassini and M. tredecassini), and those in the whole Decula (M. septendecula and M. tredecula). The whole mitochondrial sequence (1,832 bp) was treated as one locus. For each group, we attempted 7–12 runs to find the appropriate priors of t (time of population splitting), q (effective population size), and m (migration rate), with 104 or 105 burn-in generations and 104 to 106 post–burn-in generations. The mutation rate of the mitochondrial locus was set to 2.9 × 10−7 per locus per year based on the substitution rate 0.016 per site/my estimated in a BEAST analysis with the entire mitochondrial sequence data (see below).
Demographic History.
To explore the demographic history of each species group, we applied the BSP analysis (43) with the mitochondrial COI-II sequence data by using programs BEAST version 1.6.2 and TRACER version 1.5. Sequence datasets were prepared for Decim (n = 329), Cassini (n = 238), and Decula (n = 165). The whole mitochondrial sequence (1,832 bp) was treated as one locus. For each dataset, the BEAST run was conducted for 108 generations sampling every 5,000th generation. Substitution models selected by MrModeltest version 2.3 were HKY+I for Decim and Cassini datasets, and HKY for Decula dataset. A strict clock model was used with a clock rate of 0.016, which was estimated in a BEAST analysis for the entire mitochondrial dataset also with a strict clock model and a Magicicada MRCA node age prior of the normal distribution with mean 4.16 Mya and SD = 0.3943 Mya (108 generations with sampling every 104-th generation; 106 burn-in generations). After the BEAST skyline plot analysis, effective population size through time was reconstructed by using TRACER version 1.5, discarding initial 107 generations as burn-in. We also conducted an extended BSP analysis (44) by using BEAST version 1.7.2 using Cal, EF1-a, and Wg for Decim; EF1-a and Wg for Cassini; and Cal and EF1-a for Decula, in addition to COI-II data, excluding nuclear genes that showed no intraspecies group sequence difference. The number of sequences for each species group was 25, 18, and 19 for Decim, Cassini, and Decula, respectively. Because sequence difference in each nuclear gene was fairly small, nuclear gene trees were linked, and in Cassini, substitution and clock models were linked between the two nuclear genes. Substitution models selected by MrModeltest version 2.3 were HKY+G for COI-II, HKY for Cal, and GTR+I for EF1-a and Wg. A strict-clock model was used for each gene partition, and the clock rate of 0.016 was applied to COI-II gene as in the previous BSP analysis. For each species group, the BEAST run was conducted for 50 million generations sampling every 5,000th generation. Effective population size through time was reconstructed by using TRACER version 1.5, discarding initial 5 million generations as burn-in.
Data Availability
Data deposition: The sequences reported in this paper have been deposited in the DNA Data Bank of Japan (accession nos. AB740543–AB740917).
Acknowledgments
We thank David Marshall for help with field work and comments on phylogenetic methods. Many other people helped collect Magicicada specimens over the 30-y period and have been acknowledged in previous C.S. laboratory publications. This study was supported by Japan Society for Promotion of Science Grants-in-Aid 22255004 and 22370010 (to J.Y.), and 23128507 (to T.S.); Ministry of Education, Culture, Sports, Science, and Technology (the Global Center of Excellence program A06), Japan; and National Science Foundation Division of Environmental Biology Grants 8411083, 8509164, 9196213, 9812779, 9807113, 9982039, 0720664, and 0955849 (to C.S.).
Supporting Information
Supporting Information (PDF)
Correction to the Supporting Information
- Download
- 1.28 MB
st01.doc
- Download
- 351.00 KB
st02.doc
- Download
- 63.00 KB
st03.doc
- Download
- 60.00 KB
st04.doc
- Download
- 49.50 KB
References
1
CL Marlatt, The periodical cicada. USDA Bureau Entomol Bull 71, 1–183 (1907).
2
RD Alexander, TE Moore, The evolutionary relationships of 17-year and 13-year cicadas, and three new species (Homoptera, Cicadidae, Magicicada). Misc Pub Mus Zool Univ Michigan 121, 1–59 (1962).
3
M Lloyd, HS Dybas, The periodical cicada problem I. Population ecology. Evolution 20, 133–149 (1966).
4
M Lloyd, HS Dybas, The periodical cicada problem. II. Evolution. Evolution 20, 466–505 (1966).
5
RM May, Periodical cicadas. Nature 277, 347–349 (1979).
6
C Simon, Evolution of 13- and 17-year periodical cicadas (Homoptera: Cicadidae: Magicicada). Bull Entomol Soc Am 34, 163–176 (1988).
7
K Heliövaara, R Väisänen, C Simon, Evolutionary ecology of periodical insects. Trends Ecol Evol 9, 475–480 (1994).
8
KS Williams, C Simon, The ecology, behavior, and evolution of periodical cicadas. Annu Rev Entomol 40, 269–295 (1995).
9
J Yoshimura, The evolutionary origins of periodical cicadas during ice ages. Am Nat 149, 112–124 (1997).
10
PR Grant, The priming of periodical cicada life cycles. Trends Ecol Evol 20, 169–174 (2005).
11
RT Cox, CE Carton, Pleoclimatic influences in the evolution of periodical cicadas (Insecta: Homoptera: Cicadidae: Magicicada spp.). Am Midl Nat 120, 183–193 (1988).
12
J Yoshimura, T Hayashi, Y Tanaka, K Tainaka, C Simon, Selection for prime-number intervals in a numerical model of periodical cicada evolution. Evolution 63, 288–294 (2009).
13
Y Tanaka, J Yoshimura, C Simon, JR Cooley, K Tainaka, Allee effect in the selection for prime-numbered cycles in periodical cicadas. Proc Natl Acad Sci USA 106, 8975–8979 (2009).
14
C Simon, Evolution of periodical cicadas: Phylogenetic inferences based on allozymic data. Syst Zool 28, 22–39 (1979).
15
AP Martin, C Simon, Anomalous distribution of mitochondrial and nuclear gene markers in periodical cicadas. Nature 336, 247–249 (1988).
16
A Martin, C Simon, Differing levels of among-population divergence in the mitochondrial DNA of periodical cicadas related to historical biogeography. Evolution 44, 1066–1080 (1990).
17
C Simon, et al., Genetic evidence for assortative mating between 13-year cicadas and sympatric “17-year cicadas with 13-year life cycles” provides support for allochronic speciation. Evolution 54, 1326–1336 (2000).
18
JR Cooley, C Simon, DC Marshall, K Slon, C Ehrhardt, Allochronic speciation, secondary contact, and reproductive character displacement in periodical cicadas (Hemiptera: Magicicada spp.): Genetic, morphological, and behavioural evidence. Mol Ecol 10, 661–671 (2001).
19
JR Petit, et al., Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. Nature 399, 429–436 (1999).
20
DE Soltis, AB Morris, JS McLachlan, PS Manos, PS Soltis, Comparative phylogeography of unglaciated eastern North America. Mol Ecol 15, 4261–4293 (2006).
21
DC Marshall, JR Cooley, Reproductive character displacement and speciation in periodical cicadas, with description of new species, 13-year Magicicada neotredecem. Evolution 54, 1313–1325 (2000).
22
DC Marshall, JR Cooley, BR Hill, Developmental plasticity of life-cycle length in thirteen-year periodical cicadas (Hemiptera: Cicadidae). Ann Entomol Soc Am 104, 443–450 (2011).
23
RT Cox, CE Carton, Evidence for genetic dominance of the 13-year life cycle in periodical cicadas (Homoptera: Cicadidae: Maigicicada spp.). Am Midl Nat 125, 63–74 (1991).
24
RT Cox, CE Carlton, A comment on gene introgression versus en masse cycle switching in the evolution of 13-year and 17-year life cycles in periodical cicadas. Evolution 57, 428–432, discussion 433–437 (2003).
25
Y Nariai, et al., Life cycle replacement by gene introduction under an allee effect in periodical cicadas. PLoS ONE 6, e18347 (2011).
26
MJ West-Eberhard Developmental Plasticity and Evolution (Oxford Univ Press, New York, 2003).
27
MS Moulds, An appraisal of the higher classification of cicadas (Hemiptera: Cicadoidea) with special reference to the Australian fauna. Rec Aust Mus 57, 375–446 (2005).
28
R Lanfear, B Calcott, SYW Ho, S Guindon, Partitionfinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol 29, 1695–1701 (2012).
29
G Jobb, A von Haeseler, K Strimmer, TREEFINDER: A powerful graphical analysis environment for molecular phylogenetics. BMC Evol Biol 4, 18 (2004).
30
M Clement, D Posada, KA Crandall, TCS: A computer program to estimate gene genealogies. Mol Ecol 9, 1657–1659 (2000).
31
P Casgrain, P Legendre The R-package for multivariate and spatial analysis v. 4.0 5d-user’s manual (Univ de Montreal, Montreal, 2004).
32
L Excoffier, LG Laval, S Schneider, ARLEQUIN v. 3.0: An integrated software package for population genetics data analysis. Evol Bioinform Online 1, 47–50 (2005).
33
AJ Drummond, A Rambaut, BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7, 214 (2007).
34
A Papadopoulou, I Anastasiou, AP Vogler, Revisiting the insect mitochondrial molecular clock: The mid-Aegean trench calibration. Mol Biol Evol 27, 1659–1672 (2010).
35
SYW Ho, MJ Phillips, A Cooper, AJ Drummond, Time dependency of molecular rate estimates and systematic overestimation of recent divergence times. Mol Biol Evol 22, 1561–1568 (2005).
36
JAA Nylander MrModeltest v2 (Uppsala Univ, Uppsala, 2004).
37
P Vos, et al., AFLP: A new technique for DNA fingerprinting. Nucleic Acids Res 23, 4407–4414 (1995).
38
X Vekemans AFLP-SURV version 1.0 (Univ Libre de Bruxelles, Brussels, 2002).
39
Y Benjamini, Y Hochberg, Controlling the false discovery rate: A practical and powerful approach to multiple testing. J Roy Stat Soc B (Methodol) 57, 289–300 (1995).
40
J Felsenstein, PHYLIP - phylogeny inference package (Version 3.2). Cladistics 5, 164–166 (1989).
41
R Nielsen, J Wakeley, Distinguishing migration from isolation: A Markov chain Monte Carlo approach. Genetics 158, 885–896 (2001).
42
J Hey Documentation for IMa2 (Rutgers Uni, New Brunswick, NJ, 2011).
43
AJ Drummond, A Rambaut, B Shapiro, OG Pybus, Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol 22, 1185–1192 (2005).
44
J Heled, AJ Drummond, Bayesian inference of population size history from multiple loci. BMC Evol Biol 8, 289 (2008).
Information & Authors
Information
Published in
Classifications
Copyright
Freely available online through the PNAS open access option.
Data Availability
Data deposition: The sequences reported in this paper have been deposited in the DNA Data Bank of Japan (accession nos. AB740543–AB740917).
Submission history
Published online: March 18, 2013
Published in issue: April 23, 2013
Keywords
Acknowledgments
We thank David Marshall for help with field work and comments on phylogenetic methods. Many other people helped collect Magicicada specimens over the 30-y period and have been acknowledged in previous C.S. laboratory publications. This study was supported by Japan Society for Promotion of Science Grants-in-Aid 22255004 and 22370010 (to J.Y.), and 23128507 (to T.S.); Ministry of Education, Culture, Sports, Science, and Technology (the Global Center of Excellence program A06), Japan; and National Science Foundation Division of Environmental Biology Grants 8411083, 8509164, 9196213, 9812779, 9807113, 9982039, 0720664, and 0955849 (to C.S.).
Notes
This article is a PNAS Direct Submission.
See Commentary on page 6620.
Authors
Competing Interests
The authors declare no conflict of interest.
Metrics & Citations
Metrics
Citation statements
Altmetrics
Citations
Cite this article
110 (17) 6919-6924,
Export the article citation data by selecting a format from the list below and clicking Export.
Cited by
Loading...
View Options
View options
PDF format
Download this article as a PDF file
DOWNLOAD PDFLogin options
Check if you have access through your login credentials or your institution to get full access on this article.
Personal login Institutional LoginRecommend to a librarian
Recommend PNAS to a LibrarianPurchase options
Purchase this article to access the full text.