Introduction

Periodical cicadas (Magicicada spp.) of the eastern United States are remarkable for their incredibly long, prime-numbered life cycles, dense populations and almost perfectly synchronized mass emergences 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15. Meticulous surveys of "broods," populations composed of multiple species that emerge in a single year, have been carried out, and these surveys have revealed that the distribution of broods is largely parapatric 1,16,17,18. Different year classes are geographically isolated or at most coexist in the peripheries of their boundaries 1,11,16,19,20,21. Studies also revealed a few cicadas emerging off-schedule, commonly called stragglers, in both 13- and 17-year broods 22. Stragglers of 17-year broods commonly appear four years ahead of schedule, whereas four-year delays are observed in 13-year broods 11,22,23. However, field studies have found stragglers occurrence at other intervals for both 17- and 13-year broods 23. While the definite cause of straggling is not yet confirmed 11,22, studies have suggested weather-related triggers, malnutrition, nymphal overcrowding and genetic mutations caused by 13- and 17-year hybridization as possible causes 2,4,5,20,21,24,25,26. However, due to the “predator foolhardy” behavior of periodical cicadas, afforded by their massive population, most stragglers do not survive due to heavy predation and other causes 9,10,20,21,27.

Stragglers are believed to play an important part in brood formation (see Fig. 6E of 10, which was most recently updated in Fig. 8 of 18). Indeed, if population densities are high enough for predator satiation, stragglers can potentially create another temporally isolated population 10,20,21,22. However, despite the common occurrence of stragglers, the recent confirmation of nonoverlapping distributions, which resemble “puzzle pieces,” suggests some mechanism or competitive interaction that limits brood overlaps 1,18. This mechanism limits the stragglers from establishing a new self-sustaining population within the territory of the main brood. If stragglers cannot invade the territory of the main brood, then the geographic coexistence of two broods (in the same area) becomes impossible.

In the past, many mathematical studies have attempted to explain the remarkable traits of periodical cicadas. Periodicity was shown to have developed during the ice ages due to extending the nymphal stage (nymphs cease to grow if temperature is below a threshold) and Allee effects 7,28. Hybridization was suggested to lead to reduced population densities since portions of the population would emerge at different periods and would have difficulty surviving due to heavy predation and other causes. Then, the prime-numbered life cycle would have been selected for due to rare hybridization instances with other life cycles 2,13. Synchronized emergence for a life cycle longer than 10 years was also shown to be a product of predation and the limited carrying capacity of the environment 29. Recently, the nonoverlapping distribution of broods has been explained by mathematical models assuming predation on adults and competition among nymphs of the same and different broods 30,31. By assuming either a steady annual rate of straggler emergence or a single massive straggler emergence and different forms of competitive interactions, they determined the conditions that lead to the replacement of the original brood. However, the plausibility of direct interbrood competition remains uncertain since the nymphs, which are root xylem feeders 6,32, move deeper underground as they grow and avoids interference 20,33,34. Furthermore, there is little empirical evidence to ascertain competition between nymphs of different ages and to determine the functional form of competition 30.

We hypothesize that predation on a single year, together with a limited carrying capacity of the environment, would provide the segregation mechanism for maintaining the nonoverlapping distribution of different broods, despite the occurrence of stragglers, without assuming inter-brood competition. Our model is built with reference to the predation model of Bulmer for \(k\)-periodical insects 14, which is given as follows:

$$\begin{array}{ll}{n}_{j+1}\left(t+1\right)={F}^\frac{1}{k}{n}_{j}(t)\text{exp}\left(-{\sum }_{i}{\beta }_{ij}{n}_{i}\left(t\right)\right),& 0\le j\le k-2\\ {n}_{0}\left(t+1\right)={F}^\frac{1}{k}{n}_{k-1}\left(t\right)\text{exp}\left(-\sum {\beta }_{i,k-1}{n}_{i}\left(t\right)\right)\left(1-\text{exp}\left(-\frac{{n}_{k-1}\left(t\right)}{c}\right)\right),& \text{otherwise}\end{array},$$
(1)

where \({n}_{j}(t)\) is the population of the insect of age \(j\) (\(j=\text{0,1},\dots k-1\)) in year \(t\), \(F\) is the growth rate in the absence of competition, \({\beta }_{ij}\) is the competition coefficient between year-classes with ages \(i\) and \(j\), and \(c\) is the intensity of predation at age \(k-1\) (we inverted \(c\) for convenience sake). This model assumes predation on a single year and a Holling Type II functional response of predators, that is, the portion of preys eaten when preys are abundant is less than when preys are scarce due to predator satiation. This is represented by the term \(\text{exp}\left(-\frac{{n}_{k-1}\left(t\right)}{c}\right).\) The term \(\text{exp}\left(-{\sum }_{i}{\beta }_{ij}{n}_{i}\left(t\right)\right)\) is a density-dependent term, which incorporates competition within and between year-classes, and prevents explosion of the insect population. If there is no competition between (within) year-classes, then \({\beta }_{ij}=0\) (\({\beta }_{ii}=0\)). The model has unstable and stable equilibria if the predation intensity \(c\) is not too strong. Depending on whether the initial population is greater or less than the unstable equilibria, the population approaches the stable equilibria or goes extinct. Thus, the predation intensity plays a crucial role in this model. Bulmer concluded that the importance of predation is in reinforcing the tendency of competition between year-classes to lead to periodical behavior, but that it will not cause periodical behavior in the absence of such competition 14.

In this study, we maintain the predation term and within year-class competition and eliminate competition between year-classes. We demonstrate that predation provides the segregation mechanism by showing the extinction of an additional insect population (stragglers population). The proposed mechanism arises from (i) a reduction in the difference between stable and unstable equilibria as predation intensity approaches the critical value beyond which the population becomes extinct, and (ii) population imbalance between the main and straggler broods.

Model and results

Let \({N}_{0}\) be the population size of the adult cicadas. The population size in the next year is given by \({N}_{1}=F\left(1-\frac{{N}_{0}}{K}\right){N}_{0}\), where \(F\) is the growth rate (power of increase per generation in Bulmer 25) and \(K\) is the carrying capacity for the adult cicadas. This term imitates the competition within year-class term in Eq. (1). We assume predation on the first instar since most nymphal mortality is experienced at the soil surface during the first two years 25. Furthermore, predator satiation is also likely to occur due to their massive population. In accordance with Bulmer’s model (Eq. (1)), we describe the predation effect by means of a parameter \(c\), the predation intensity. The population size of the second instar \({N}_{2}\) is given by \({N}_{2}=\left(1-\text{exp}\left(-\frac{{N}_{1}}{c}\right)\right){N}_{1}\). The results are not affected essentially if predation occurs at any other age (see Supplementary Materials 1). The risk of mortality is little to none afterward and the cicadas maintain their population until emergence 25. To summarize, the model is given by

$${N}_{t+1}=\left\{\begin{array}{ll}F\left(1-\frac{{N}_{t}}{K}\right){N}_{t}& \text{mod}\left(t,17\right)=0\\ \left(1-\text{exp}\left(-\frac{{N}_{t}}{c}\right)\right){N}_{t}& \text{mod}\left(t,17\right)=1\\ {N}_{t}& \text{otherwise}\end{array}\right.,$$
(2)

where \(\text{mod}(t,17)\) is the remainder of \(t\) divided by \(17\). This model has a trivial solution \({N}_{t}=0\), i.e., extinction. If the intensity \(c\) is not too strong, there are stable and unstable equilibrium densities, which are obtained from the equilibrium condition \({N}_{17}={N}_{0}\) (see Supplementary Materials 2). The equilibrium densities \({N}_{0}({=N}_{17})\) are plotted against \(c\) in Fig. 1 (\(F=3, K=400\)). If the initial population size is above (below) the unstable equilibrium, the population approaches the stable equilibrium (goes extinct). If the predation intensity is beyond a critical value (\(c>390\)), the equilibria cease to exist, and the population will always become extinct. As the intensity \(c\) approaches the threshold from below, stable and unstable equilibria merge to the critical point. Thus, the main brood at the stable equilibrium is more likely protected from stragglers if the intensity is close to the critical value. On the contrary, it is comparatively easier for stragglers to survive and reach an equilibrium if the intensity is weak.

Fig. 1
figure 1

Stable and unstable equilibria plotted against predation intensity \((F=3, K=400)\). If the initial population size is above (below) the unstable equilibrium, the population approaches the stable equilibrium (goes extinct). If the predation intensity is beyond a critical value (\(c>390\)), the equilibria cease to exist, and the population will always become extinct.

In Fig. 2, we present the results of numerical simulation for \(c=250\) and 350, where the initial population is \({N}_{0}=100\) and 50 (\(F=3, K=400)\). Being the only case where the initial population size is below the unstable equilibrium (as shown Fig. 1), extinction occurs only in the case of \(c=350\) and \({N}_{0}=50\) (Fig. 2d), whereas stable equilibria are reached after \(t=187\) in the other three cases (Fig. 2a–c). As the predation intensity \(c\) increases from 250 (Fig. 2a,c) to 350 (Fig. 2b), the portion of captured preys \(\text{exp}\left(-\frac{{N}_{17n+1}}{c}\right)\) increases from 0.3 to 0.4, respectively. Thus, the adult population in equilibrium decreases from 209 to 163, while the amplitude of oscillation (\({N}_{17n+2}-{N}_{17n+1}\)) increases. Supplementary Materials 3 provides illustrations for 13-year broods.

Fig. 2
figure 2

Simulated population dynamics of a single brood. (a) \(c=250\) and \({N}_{0}=100\). (b) \(c=350\) and \({N}_{0}=100\). (c) \(c=250\) and \({N}_{0}=50\). (d) \(c=350\) and \({N}_{0}=50\) (\(F=3, K=400\)). Being the only case where the initial population size is below the unstable equilibrium, extinction occurs only in the case of \(c=350\) and \({N}_{0}=50\), whereas stable equilibria are reached after \(t=187\) in the other three cases.

After the population of the main brood \({N}_{t}\) reaches a stable population at time \(T>0\), we introduce a straggler brood \({J}_{t}\). We use a model similar to Eq. (2) to describe the population of the stragglers that emerges \(y\) years after the main brood, that is,

$${J}_{t+1}=\left\{\begin{array}{ll}F\left(1-\frac{{J}_{t}}{K}\right){J}_{t}& mod\left(t,17\right)=y\\ \left(1-\text{exp}\left(-\frac{{J}_{t}}{c}\right)\right){J}_{t}& mod\left(t,17\right)=y+1\\ {J}_{t}& \text{otherwise}\end{array}, t>T\right..$$
(3)

In Fig. 3, adult stragglers \({J}_{T+y}=30\), 40, 100 and 170 are introduced \((c=250, F=3, K=400, T=187, y=13)\). In so far as the straggler population is small, i.e., below the threshold (unstable equilibrium in Fig. 1), they did not survive and are eliminated after some generations (Fig. 3a,b). However, given a sufficient initial population, the stragglers can reach a stable equilibrium (Fig. 3c,d). The two broods both reach stable populations owing to the lack of direct competition in the present model (Fig. 3c). In the case of a massive off-scheduled emergence shown in Fig. 3d, the main brood population would only be left with around 40 adults, which results in its extinction.

Fig. 3
figure 3

Simulated population dynamics of main and straggler broods. (a) The initial population size of a straggler brood is \({J}_{T+y}=60\). (b) \({J}_{T+y}=\) 70. (c) \({J}_{T+y}=\) 80. (d) \({J}_{T+y}=\) 90 \((c=250, F=3, K=400, T=187, y=13)\). In so far as the straggler population is small, they did not survive and are eliminated after some generations. However, given a sufficient initial population, the stragglers can reach a stable equilibrium.

Discussion

Many studies have tried to explain how periodical cicadas developed their extraordinary characteristics. For instance, it was suggested that common temperate cicadas evolved to become the periodical cicadas 28. Suggestions on the possible cause of the long life cycle include the delayed development of nymphs during the Pleistocene era, predator escape, and selection for large body size 7,10,20,25,28. Periodicity and the fixed 17- and 13-year life cycles would have been selected for due to rare hybridization instances with other life cycles and Allee effect 2,13,15. Hybridization would have resulted in reduced population densities since portions of the population would emerge at different periods, and these emergences would have been eliminated by predators or other causes. Finally, the remarkable abundance was suggested to be an effect of periodicity, synchrony, and long life cycle 9,10.

This study provides an explanation for the nonoverlapping distribution of periodical cicada broods through mathematical models and simulations. Periodical cicadas are known for “predator-foolhardiness” and do not defend themselves against predators 10. However, the massive population of these insects allows for a still sizeable number to escape predation. Our results have shown that predation can provide a mechanism for maintaining the nonoverlapping distribution of broods despite the common occurrence of stragglers. Under suitable circumstances and predation intensity, the persistence of a brood depends on the initial population size of the stragglers, but it does not depend on the timing of their introduction. If the initial population of the stragglers is above the threshold (solid line in Fig. 1), the stragglers can reach a stable equilibrium and coexist with the main brood. However, in case of massive stragglers emergence when the population of the main brood falls below the threshold, the main brood is eliminated, and brood replacement occurs. Under strong predation intensity, both the main brood and stragglers can be eliminated. In theory, these results should generate testable predictions concerning the temporal separation or overlap of broods’ border; however, to date, such predictions have been difficult to test empirically.

In 17-year periodical cicadas, past studies have revealed that four-year accelerated stragglers are most common, while one-year accelerated ones are also observed, especially in the Midwest 11,22,23. However, stragglers of the other year classes (e.g., delayed and two- or three-year accelerated off-scheduled emergences) have been rarely reported. There are two competing hypotheses for the four-year accelerated stragglers. The first hypothesis is genetic mixing between 17- and 13-year cicadas 26,35,36. Sympatric broods of homozygous 13- and 17-year cicadas emerge together and crossbreed every 221 years. The F1 hybrids are heterozygous and emerge together with the homozygous broods, while some emerge in other years (both early and late) 2,26. The F2 generation will include homozygous recessive individuals. These homozygous recessive individuals may produce a new brood, which emerges four years earlier or later than the parent 17-year brood. While it is yet to be established, the longer 17-year gene is likely to be recessive. According to the other hypothesis, the four-year acceleration occurs due to some intrinsic or extrinsic causes other than hybridization 22. Under this mechanism, four-year stragglers emerge every generation spontaneously, while these stragglers may leave no offspring.

There are suggestions that 17-year cicadas emerging four years before the main brood, with a dense enough population to satiate the predators, can establish a new brood 1,20,21,22,24. This suggestion was birthed from observations that broods that are four years apart coexist or live in the peripheries of each other 1,16,20,24. For instance, Broods XIX and XXIII in Missouri 19, and Broods XIV, X, IX, V and I in Long Island, New York 20,24 are suggested to at least share boundaries. Our model has shown that a single main brood can birth an independent self-sustaining brood given that the stragglers reach the threshold presented in Fig. 1. Alternatively, stragglers can join the (on-time) emergence of a neighboring brood and be afforded the numbers to escape predation. In either case, if the main brood emerges later with sufficient numbers to satiate predators, both broods are expected to reach stable populations and coexist (Fig. 3c). However, in case the remaining population of the main brood drops below the threshold due to massive straggler emergence, brood replacement would occur (Fig. 3d). It takes several generations or centuries before both broods reach stable populations or before the main brood population is completely eradicated, which would make empirical confirmation difficult. For the time being, it would seem that two broods are coexisting, albeit one brood with a higher population. The current model could not explain the bias of four-year differences in neighboring broods.

Birds are the most observed predators of adult cicadas 8,37. Avian predation is hypothesized to be a major driver of the synchronous emergence of periodical cicadas and has been observed to eliminate adult cicadas that emerge out of schedule 8. However, other inconspicuous ground predators have also been observed to prey on emerging fifth instar nymphs, such as spiders and wasps 38. Whereas nymphs are too small for birds to prey on, such insects can play a part in eliminating stragglers and maintaining the nonoverlapping distributions. Ants have also been observed to prey on freshly hatched nymphs that fell on the ground 9,25,39 (see also Supplementary Materials 4). Ants are social insects and orient their nest mates to the ___location of food sources 40,41. Ants commonly exhibit colonial hunting when the food resource (cicada nymphs) is abundant, i.e., many workers come out for foraging 42,43. Ant colonies are the most abundant in forest environments among all predators 44. A previous study has estimated that over 90% of periodical cicada mortality is experienced during the hatching and the first instar nymph stage (while digging into the soil) and almost no mortality thereafter 25. This disqualifies moles as the main source of mortality since first instar nymphs are too small for moles to prey on 9. We suggest the possibility of predation by ant colonies as the source of this mortality in periodical cicadas. It is left for a future study to confirm if ants are the primary first instar predators and to identify other possible predators.

Our findings offer intriguing insights into Lloyd and Dybas’ 10 discussion of the evolution of periodicity, and how different sources of mortality lead to different patterns regarding cicadas coming out before the brood emergence (“leading” emergers) and cicadas coming out after the brood emergence (“lagging” emergers). Lloyd and Dybas hypothesized that “leading” off-cycle emergences would generally be favored except in the case of a specialized, synchronized parasitoid or predator, which would impose strong selection against such early emerging cicadas. Such emergence-timing dependence is not observed in the present model, for the predation intensity is assumed not to vary in response to the population dynamics. However, it is possible that predation intensity varies somehow in response to mass emergence. If the intensity increases for years after mass emergence, lagging emergence should be eliminated more effectively than leading emergence.

Direct competition between nymphs has been considered a viable mechanism 14,30,31. However, empirical studies show that no competition between nymphs of different sizes (brood) takes place among coexisting three species 45. White and Lloyd have found that nymphs move deeper underground as they grow 34. Consequently, competition for feeding sites becomes more unlikely. Furthermore, three species coexisting in a single brood do not compete with each other 46. As a matter of fact, the present study was motivated from this observation of the unlikelihood of competition.

The evolution of periodical cicadas is a very unique event in the ice ages, consisting of two sequential processes 2,28. All common cicadas in temperate regions have size-dependent maturity, where life cycles vary depending on available food (water) resources 47. The genus Magicicada has the longest known fixed life cycles, significantly longer than the typical life cycles of temperate cicadas (at most ten years) 9,47. The elongation of life cycles and resulting periodicity should be caused by geological cooling during the ice ages in the Pleistocene 7,28. The selection of 13- and 17-year (prime numbered) life cycles should have occurred afterward, either during the inter-glacial or post-glacial periods 13,15,28. Molecular phylogeny of periodical cicadas demonstrates that the present Magicicadas were established in the Pleistocene 12.

The confirmation of nonoverlapping distribution of broods is relatively recent than the other characterizing features of periodical cicadas 1,16,18 and with this comes the investigation of its probable causes 30,31. This feature is perhaps not yet appreciated in the original predation model 14 but the results are similar in a sense that the equilibrium depends on the intensity of predation of the initial size of the population. There is still much to be understood regarding how periodical cicadas achieve nonoverlapping distribution of broods, and there are even more questions on the mechanisms driving stragglers and off-schedule emergences. Future studies can be done to generalize the predation intensity \(c\) or to consider other models that can account for the numerical response of predators. The current model also only allows a single introduction of stragglers, which follows the life cycle of the main brood. Future studies can focus on multiple introductions of stragglers and/or contraction or elongation of life cycles. Such models may better explain the four-year differences between neighboring broods and better support brood replacement theories 24. The bias in four-year and one-year acceleration of 17-year cicadas and four-year delay of 13-year cicadas still needs to be studied.

Methods

We employ difference equations to examine if predation can provide a separation mechanism for broods of 17-year periodical cicadas. To this effect, we modify Eq. (1) to reflect the life cycle of the main brood (Eq. 2) as follows:

  1. i.

    Adult cicadas emerge and reproduce every 17 years, and die thereafter. If \(\text{mod}\left(t,17\right)\) denotes the remainder when \(t\) (in years) is divided by 17, then, the main brood emerges and reproduces when \(\text{mod}\left(t,17\right)=0\). For simplicity, we assume that the cicadas spent a year above-ground, from emergence to hatching.

  2. ii.

    We incorporate a carrying capacity of the environment for the adult cicadas. Periodical cicadas are known for their incredible abundance, which can lead to competition for food, mates, and suitable twigs for ovipositing 9,10,20. This replaces the within-brood competition term in Eq. (1).

  3. iii.

    We assume that predation occurs during the first instar only, \(\text{mod}\left(t,17\right)=1\), when the highest mortality of nymphs is observed 25,27. We maintain the Holling Type II functional response, i.e., we assume predator satiation.

  4. iv.

    Since there are little to no mortality risks after the first instar 25, we assume that the cicadas maintain their population until their next emergence.

  5. v.

    We remove the between-brood competition term since nymphs are found to move deeper underground as they grow and avoid interference 20,33,34.

Numerical simulations are conducted to parameterize Eq. (2) (codes are written using Python 3.9, see Data Availability Statement). We only consider Eq. (2) in parameterizing the model since the straggler brood population is introduced only after the required equilibrium for the main brood population is reached. In case the formula return a negative value for \(N\) or \(J\), we record the population as 0. We test the model for \(F=\text{1,2},\text{3,4},5\), \(K=100, 200, 300, 400, 500\) and \(c\in [\text{0,500}]\) for 40 generations. The values for the parameters are chosen so that there are 200–300 adult cicadas at equilibrium 3,25,27. This (periodical) equilibrium oscillates from a minimum population of 200–300 cicadas every \(\text{mod}\left(t,17\right)=2\) after predation events and the rest of the cicada life cycle to a maximum every \(\text{mod}\left(t,17\right)=1\) after the emergence-reproduction period. For illustration, Figs. 1, 2 and 3 are generated using \(F=3, K=400\).

To investigate if stragglers can disrupt an already established brood, the introduction of stragglers is done only after an equilibrium for the main brood population is reached at time \(T>0\), \(\text{mod}\left(T,17\right)=0\). The life cycle of straggler brood is the same as the main brood except the straggler brood emerges \(y\) years after the main brood, \(\text{mod}\left(t,17\right)=y\), and predation events whenever \(\text{mod}\left(t,17\right)=y+1\) (Eq. 3). In Fig. 3, we choose \(T=187\) (11 generations) and \(y=13\). The stragglers are off-scheduled emergence from the main brood. In effect, the main brood population is reduced by the number of stragglers that emerge. Stragglers are only introduced once and their offspring emerge every 17 years thereafter if they survive. Periodical cicada densities have no upper or lower limits below which they cease reproducing 9,10,25,27. Thus, nymphs that survived above-ground predation would reproduce and have the potential to establish a new self-sustaining brood. We observe the dynamics of the broods 20 generations after the stragglers are introduced.