A second view on the evolution of flight in stick and leaf insects (Phasmatodea)

Background The re-evolution of complex characters is generally considered impossible, yet, studies of recent years have provided several examples of phenotypic reversals shown to violate Dollo’s law. Along these lines, the regain of wings in stick and leaf insects (Phasmatodea) was hypothesised to have occurred several times independently after an ancestral loss, a scenario controversially discussed among evolutionary biologists due to overestimation of the potential for trait reacquisition as well as to the lack of taxonomic data. Results We revisited the recovery of wings by reconstructing a phylogeny based on a comprehensive taxon sample of over 500 representative phasmatodean species to infer the evolutionary history of wings. We additionally explored the presence of ocelli, the photoreceptive organs used for flight stabilisation in winged insects, which might provide further information for interpreting flight evolution. Our findings support an ancestral loss of wings and that the ancestors of most major lineages were wingless. While the evolution of ocelli was estimated to be dependent on the presence of (fully-developed) wings, ocelli are nevertheless absent in the majority of all examined winged species and only appear in the members of few subordinate clades, albeit winged and volant taxa are found in every euphasmatodean lineage. Conclusion In this study, we explored the evolutionary history of wings in Phasmatodea and demonstrate that the disjunct distribution of ocelli substantiates the hypothesis on their regain and thus on trait reacquisition in general. Evidence from the fossil record as well as future studies focussing on the underlying genetic mechanisms are needed to validate our findings and to further assess the evolutionary process of phenotypic reversals. Supplementary Information The online version contains supplementary material available at 10.1186/s12862-022-02018-5.


Introduction 26
Active flight is considered the key innovation to have driven lineage diversification in animals 27 and has allowed insects to become the most species-rich group on Earth 1-3 . The evolution of the 28 insect wing ~400 million years ago (mya) 2,4 and the associated enhanced abilities to disperse 29 and access unreachable resources are directly linked to the remarkable radiation and success of 30 winged insects (Pterygota). Moreover, wings have undergone various modifications to 31 additionally or alternatively serve other functions such as thermoregulation, mate choice and 32 courtship, crypsis and defensive strategies 1,5 . However, despite the numerous advantages, wings 33 have been repeatedly reduced across all pterygote groups with several lineages being 34 completely wingless such as heel-walkers (Mantophasmatodea), lice (Phthiraptera) or fleas 35 (Siphonaptera) 1,6-8 . Paradoxically, also the loss of flight is proposed to be in direct correlation 36 with increased diversification rates and thus also a driver of speciation [9][10][11][12] . 37 Flightlessness appears to occur when the selection for aerial dispersal is relaxed as in habitats 38 with environmental stability or with unfavourable conditions for flight such as high winds or 39 cold temperatures 6,13,14 . As a consequence, flightless taxa may show increased fecundity because 40 resources previously invested in maintaining the energetically costly flight apparatus can now 41 be allocated to reproduction 6,15,16 . The underlying trade-off between dispersal and reproduction 42 has been repeatedly demonstrated for females 17-21 and males 22-24 , albeit the loss of flight is 43 generally more common in females, often resulting in wing dimorphic species with volant 44 males 1,6,25 . In some species, the dilemma of dispersal capability versus fecundity is solved by 45 wing polyphenism, where an either winged or flightless phenotype is adopted in response to 46 specific environmental triggers 1,16,18,26-29 . Yet, flightlessness is not tantamount to the complete 47 loss of wings (aptery): A phenotype or species may exhibit shortened wing length or retain fully-48 sized wings but with reduced flight musculature 26,27 . 49 Within winged insects, the plant-mimicking lineage of stick and leaf insects (Phasmatodea) has 50 been recognised as an expedient study system to investigate the evolution of flight due to their 51 high diversity in wing states and sizes varying among closely related species and between 52 sexes 30,31 (Figure 1). Generally, forewings are modified to abbreviated and sclerotised tegmina 53 or wing pads, while the hindwings are membranous and folded neatly against the elongated and 54 often slender body at rest 2 . However, most of these herbivorous insects are flightless 1,30 , being 55 either completely wingless (apterous; Figures 1A and 1B) or short-winged 56 (micropterous/brachypterous; Figures 1C-1F), while long-winged forms (macropterous; 57 Figures 1G and 1H) may or may not be capable of ascending flight. Wings that do not sustain 58 powered flight may serve for gliding or other derived utilities such as defensive stridulation or 59 startle displays 32-34 . The latter is common in both short-and long-winged species whose wing 60 undersides or membranes show bright warning colours in contrast to the otherwise 61 inconspicuously coloured body maintaining crypsis ( Figures 1E-1H). Sexual wing dimorphism 62 occurs throughout the phasmatodean lineage and is exclusively female-biased 30,31 ( Figure 1I) 63 with the exception of two Phasmotaenia species 35 . As a result, females tend to be larger and 64 brachypterous or apterous, whereas males are smaller and either share the female's wing state 65 or have more developed wings.    Herein, we revisit the evolution of flight in stick and leaf insects by inferring a phylogenetic 115 reconstruction of wings and ocelli based on the largest taxon sample to date with >500 species 116 covering all major groups. This extremely dense taxonomic resolution appears to be necessary 117 to satisfactorily reconstruct this highly disparate character system.  Prior to reconstructing ancestral states, we fitted different models to the binary wing dataset 154 (datawings) and compared log-likelihoods and AIC values (Supplementary file 7). The ARD model 155 was determined as the best-fit model independent of whether the root state was equally 156 probable to be winged/wingless or forced to be wingless. By contrast, the models disallowing 157 the regain of wings (IRR) performed the poorest, even when the root was forced to be winged. The evaluation of the number of transition events between states when wing-regain is permitted 173 (ARD) or wing loss is irreversible (IRR) showed that when analysing the 3-state dataset under 174 the ARD model, the number of shifts was highly increased, albeit most shifts were detected 175 between micropterous and macropterous taxa. Comparing the binary dataset under ARD and 176 IRR models and thus the transitions between the presence and absence of wings revealed more

Trait correlation 214
We tested for correlation of wings and ocelli and recovered strong evidence for their correlated 215 evolution in males and females (Supplementary file 10). Specifically, models performed best 216 when ocelli were set as the depending variable implying that the evolution of ocelli was 217 depending on the evolution of wings. However, applying different root states influenced the 218 estimation resulting in varying AIC values and Akaike model weights (AICw). In one instance, 219 when the root state for females was set to wingless + no ocelli, the evolution of ocelli and wings 220 was recovered as intradependent. The wingless root state also resulted in the lowest AICw 221 values for males (AICw = 0.5136), whereas the unconstrained, winged + no ocelli and winged + 222 ocelli root states resulted in higher values (AICw = 0.6462, 0.9434 and 0.7895, respectively). The 223 correlation of wings and ocelli was further examined by fitting models of correlated evolution 224 between ocelli and the three individual wing states. While we did not find any evidence that 225 short wings and ocelli are correlated, there was statistically significant correlation recovered 226 between ocelli and long wings as well as winglessness. The second-best shift configuration additionally shows a slight rate shift for the lineage 234 including among others the Cladomorphinae and Lanceocercata. Overall, the net diversification 235 shows a gradual and constant increase in rate over time and remains comparably low only for 236 Timematodea (Supplementary file 11). The results of our trait-dependent diversification rate 237 analysis to assess the influence of wings on the diversification process show that character-238 independent HiSSE (CID-4) models were inferred to best fit our data (Supplementary file 12). 239 The best-fit model (AICw = 0.9412), the CID-4 model with four hidden states, does not include 240 the focal character (wings) and thus assumes that changes in diversification rate are 241 independent from the presence or absence of wings.  and not only for this clade, where morphology is corroborative, but also for the more diverse 280 winged lineages of Euphasmatodea. 281 Regardless of whether wings were present or absent in the common ancestor of 282 Euphasmatodea, our results show an increase in diversification rate with the start of their 283 radiation (Supplementary file 11). This rapid radiation that largely follows the Cretaceous-284 Palaeogene (K-Pg) mass extinction (~66 mya) was previously suggested to be linked to the 285 diversification of flowering plants 53 -a co-evolutionary pattern also observed in other plant-286 associated insect groups 80-82 . Our results substantiate this hypothesis, particularly, when 287 considering the results of the trait-dependent diversification analysis (HiSSE) from which we 288 concluded that flight or flightlessness might not have been the main driver of euphasmatodean 289 diversification. Alternatively, the evolution of hard-shelled eggs 52,77,83 and the acquisition of 290 endogenous pectinase genes 84 further shaping the co-evolution with their food plants are most 291 likely to have contributed to the success of the early euphasmatodean lineages. Also the 292 recurrent opportunities for colonising new land masses (i.a., the Indo-Pacific region) appears to 293 have promoted speciation 76,85 and might explain the increased diversification rate recovered for 294 Lonchodidae, Lanceocercata and relatives within Oriophasmata (square symbol, Supplementary 295 file 11). By contrast, the slight decrease in the rate for the European Bacillinae (circle symbol) 296 might be an artefact resulting from oversampling closely related taxa, most of which are only 297 represented by one single gene, but could also correspond to the ecological shift to a temperate 298 region. Given that both the presence of wings and the lack thereof were shown to potentially 299 increase diversification in insects 1,2,6,9-12 , it cannot be excluded that either condition had an equal 300 impact on the diversification of the individual euphasmatodean lineages (e.g., for the 301 predominantly macropterous and most species-rich lineage Necrosciinae versus the speciose 302 and mostly wingless Lonchodinae 66,86,87 ). Ultimately, that no significant change in rate among 303 subgroups was detected and character-independent diversification was favoured is indicative of 304 the wing state not being subject to Dollo's law of irreversible evolution 43 . 305 Within Euphasmatodea, only few lineages were recovered as ancestrally winged (e.g., 306 Pseudophasmatidae, Phylliidae), whereas in others, the winged taxa are found rather 307 subordinated among otherwise wingless forms (e.g., Lonchodinae, Diapheromerinae). In 308 contrast to the seemingly random distribution of apterous species in ancestrally winged 309 lineages, the winged species in predominantly wingless groups are closely related, indicating 310 single origins of wing re-evolution. Discrepancies may be explained by incomplete taxon 311 sampling or wrongly recovered phylogenetic position, which might diminish in a phylogenomic 312 context or when applying an even denser taxon sampling. Conversely, the multiple instances of 313 secondarily apterous taxa found among winged lineages are the result of the higher probability 314 of wing loss in comparison to that of regain 30,39,40,45 . Generally, evolutionary transitions were 315 proposed to occur rapidly and often in Phasmatodea -not only between apterous and 316 macropterous taxa but also between wing states 31 . Previous authors argued that the probability 317 or cost ratio of wing gain was considered too high compared with that of its loss 40 . In our 318 reconstruction, the estimated number of potential regains ranges between 9 and 36 events 319 (~22) with ~42 losses (27-59) versus ~76 losses (71-82) if wing loss is considered irreversible 320 ( Figure 5), hence, permitting re-evolution appears to be more parsimonious. For instance, in 321 predominantly wingless Diapheromerinae, we find 1-3 potential wing regain events versus an 322 alternative of a total of ten losses if wing recovery was impossible. Although our results clearly 323 support the hypothesis of re-evolution, we nevertheless acknowledge that the reacquisition of a 324 complex trait such as wings must be assumed less likely than its loss, and that it is possible that 325 ultimately, independent wing loss might have occurred much more often in phasmatodeans than 326 in other insect groups. However, how can we then explain the distant disjunct distribution of 327 ocelli in Phasmatodea? 328 Insect ocelli have not been intensively studied in a phylogenetic context before, but it is well-329 known that there are winged lineages within Pterygota that lack ocelli such as earwigs 330 (Dermaptera) and webspinners (Embioptera). Although both groups have undergone a unique 331 modification of wings, their specific lifestyle led to a reduced need for aerial dispersal over time 332 the Aschiphasmatidae, sister group to the remaining Euphasmatodea (=Neophasmatodea), 341 comprise besides wingless forms several fully-winged species that are capable of active flight, 342 yet, none possess ocelli. Our results show that the presence of ocelli is restricted to five 343 subordinated and completely unrelated lineages, albeit winged and volant taxa are found more 344 frequently across Phasmatodea (Figure 4). While it is plausible to assume that ocelli might have 345 been lost in winged lineages that contemporarily lost the capability of flight (i.e., micropterous 346 lineages), the high number of macropterous and volant species lacking ocelli is perplexing. The 347 alternative scenario involving the possibility of wing regain implies that secondarily winged 348 lineages re-evolved the ocelli subsequent to or simultaneously with the recovery wings, which is 349 highly supported by our analyses. Specifically, the ancestral state of Pseudophasmatidae, 350 Palophidae and Phylliidae was recovered as possessing ocelli along with wings, whereas ocelli 351 are estimated to have occurred only in subordinated clades in Lanceocercata and Necrosciinae 352 in spite of wings and flight being recovered to have a more ancestral. The necrosciine taxa 353 represent an optimal example, since ocelli occur only in one single subclade within a highly 354 diverse macropterous lineage and are otherwise only present in two unrelated species, where 355 they must have re-appeared independently (Hemisosibia incerta, Gargantuoidea triumphalis). 356 Conversely, it might be a common trend for brachypterous species within the ocelli-bearing 357 clades to reduce ocelli (as observed in other insect groups 90-93 ), albeit there are also few 358 macropterous species lacking ocelli (e.g., Creoxylus, Eurynecroscia, Paraprisopus). Necrosciinae would undergo specific modifications to improve the performance of the ocellar 382 system to compensate for the obscured horizon 62,105 , yet, there are no habitat differences 383 between ocellate and anocellate Necrosciinae. Hence, we favour the hypothesis that ocelli are 384 not plesiomorphic with copious losses but instead were re-evolved as was also proposed for 385 lineages of the predominantly anocellate beetles 64,107 . To further corroborate this assumption, 386 future studies are needed to elucidate the organisation and development of ocelli across 387 phasmatodean taxa as well as to examine the internal morphology of the head capsule to clarify 388 the potential presence of internal structures in anocellate species that facilitate ocelli 389

emergence. 390
Under the assumption of Dollo's law of irreversible evolution, the developmental genetic 391 foundation of a complex trait is lost along with the trait itself and therefore cannot be reacquired 392 because evolving the same complex structures de novo is not considered possible. However, if 393 the underlying genetic architecture for the lost trait remains conserved, reinitiation of 394 suppressed genes might result in its recovery. Moreover, evolutionary changes in gene 395 regulatory networks are proposed to have a major influence on morphological evolution 108,109 . 396 Given that the genetic components and developmental processes of a complex character were 397 conserved, evolutionary shifts in developmental gene expression patterns might be responsible 398 for the loss as well as the regain of a trait 108 . The key elements underlying insect flight have been 399 conserved over hundreds of million years of evolution and may remain available for reactivation 400 in secondarily wingless insects for a long period of time 110 , albeit it was estimated that silenced 401 genetic pathways cannot retain function for more than 10 million years 111 . Yet, regulatory 402 elements may be reinitiatable over longer periods of time 50 , especially when the genetic 403 developmental programme is underlying pleiotropic maintenance and is accommodated in 404 similar structures or tissues still present 30,47,48,50,108 . Since the genetic basis for insect wing 405 development is retained in leg formation 112,113 , genes may have been co-opted from these 406 existing expression patterns that then elicited wing re-evolution in stick and leaf insects 30 . 407 Although a gradual process of the reversal of wing loss involving the necessary deployment of 408 multiple genes must be assumed, it has been shown that musculature and innervation associated 409 with flight can be maintained in flightless phasmatodeans 114 , while in other insects flight loss 410 may be explained by degradation of these structures 27 . Therefore, we suggest that in 411 phasmatodeans apparently a latent and hitherto not understood capacity to produce wings and 412 flight-associated traits exists (or existed) that facilitated the reappearance of wings, flight and 413 ocelli to occur multiple times independently. Our hypothesis thus incorporates that the disjunct 414 distribution of these traits did not result from plesiomorphy and numerous losses but from a 415 radiation of wingless ancestors whose descendants independently regained wings (and ocelli).  We applied two approaches to infer the phylogenetic relationships: Maximum Likelihood (ML) 458 and Bayesian Inference (BI). For the ML tree inference, we used IQ-TREE 128,129 to conduct 100 459 independent tree searches using the aforementioned constraints and the best-fit partitioning 460 scheme and substitution models. For each search, we used a random starting tree, increased the 461 number of unsuccessful iterations (-nstop 500) and set the perturbation strength to 0.2. We 462 computed the log-likelihoods and at the same time tested the tree topologies to identify the best-463 scoring tree. Nodal support was subsequently calculated using the Ultrafast Bootstrap 464 approximation (UFBoot) and the SH approximate likelihood ratio test (aLRT) 130 by generating 465 10,000 replicates. UFBoot support values were plotted on the best tree, whereas for the SH-aLRT 466 test, a new phylogeny was generated by default in IQ-TREE. 467 For the BI tree, we simultaneously conducted the phylogenetic analysis with the divergence time 468 estimation in BEAST2 v. 2.6.1 131 . We used the same partitioning scheme as for the ML analysis, 469 but substitution models were selected by bModelTest v. 1.2.1 132 implemented in BEAST. We 470 linked trees and clocks across all partitions and applied the Yule tree prior and a relaxed clock 471 with an uncorrelated lognormal distribution (UCLD) and a clock-rate of 1e -7 . We time-calibrated 472 the tree using five fossils unambiguously assigned to Phasmatodea (Supplementary file 16) and 473 chose a lognormal prior distribution with the minimum age selected as the offset and log-mean 474 and log-sd set to 1.0. In the case of the two fossils from Dominican amber (Clonistria sp. and 475 Malacomorpha sp.), the calibrated node was selected to include the sister group (option "use 476 originate"). For these fossils as well as for the Phasmatodea calibration, the log-mean was set to 477 2.0 to allow a softer maximum bound. To account for potential phylogenetic incongruences that 478 might have an impact on the reconstruction of ancestral states, we performed three independent 479 analyses based on different constraints: For the first, we used the same constraints as for the ML 480 analysis (B1), for the second we additionally constrained the Heteropterygidae as the sister 481 group to all the remaining Oriophasmata (following the results of transcriptomic studies 53,54 ) 482 (B2), and for the third, we did not set any constraints for Oriophasmata and Occidophasmata 483 and thus only for Phasmatodea, Euphasmatodea and Neophasmatodea (B3). All analyses were 484 run for 300 million generations with parameter and tree sampling every 10,000 generations. (www.phasmida.speciesfile.org, www.phasmatodea.com). While we assembled the ocelli data 494 based on their presence or absence and did not distinguish between different types of ocelli (e.g., 495 number, relative size), wing states were coded as wingless (apterous, 0), partially-winged 496 (micropterous, 1) or fully-winged (macropterous, 2). The differentiation between micropterous 497 and macropterous was made arbitrarily as we defined a morph as fully-winged when the hind 498 wings exceeded the fourth abdominal segment 77 . Since some analyses require binary data, we 499 created two additional datasets by combining the three wing states based on wingless (0) versus 500 winged (1) taxa (datawings) and based on flightless (0) versus fully-winged (1) taxa (dataflight). 501 Missing data was generally coded as absent (0). For all analyses, the fossil-calibrated BEAST tree 502 with B2 constraints was used. to the binary dataset of male wing traits (datawings). Using fitDiscrete, we tested models that 520 permitted evolutionary rates between wingless and winged to be equal in both directions 521 (equal-rates, ER), different in all directions (all-rates-different, ARD) and unidirectional from 522 winged to wingless (customised irreversible model, IRR). We repeated the tests using the 523 function fitMk and additionally tested the IRR model with the root state fixed to "winged" and 524 the ARD model with the root state fixed to "wingless". Finally, for the ancestral state 525 reconstruction, we used the 3-states datasets of wings and the binary datasets of ocelli for both 526 males and females, and coded missing data with equal probabilities for each of the possible 527 states. We performed stochastic character mapping 143 with the function make.simmap 528 implemented in phytools based on the best-fit model (ARD) and ran 300 iterations of Markov 529 chain Monte Carlo (MCMC) sampling. Using the option Q="mcmc", the rate matrix Q is sampled 530 300 times from its posterior probability distribution using MCMC and stochastic maps and node 531 states are generated conditioned on each sampled value of Q. To assess whether the topology of 532 the tree had an impact on the estimated ancestral states, we also performed the analysis using 533 the alternative trees (B1 and B3 constraints). 534 Given the 300 inferred possible evolutionary histories resulting from stochastic character 535 mapping, we calculated the average number of transition events between wings states and the 536 presence/absence of ocelli. For comparison to an evolutionary hypothesis where wing/ocelli 537 regain is not permitted, we repeated simulated stochastic maps as described above but using the 538 IRR model. To be even more conservative, we also reconstructed ancestral states using both 539 models for the binary wing dataset (datawings). 540 541

Trait correlation 542
We used Pagel's binary character correlation test 144 as implemented in phytools to assess 543 potential coevolutionary dynamics between wings and ocelli in males and females. By applying 544 the ARD substitution model and using different parameters for the argument "dep.var", we 545 estimated whether wings and ocelli evolved independently, intradependently or whether the 546 evolution of wings was dependent on the evolution of ocelli or vice versa, and determined the 547 best model to explain the correlated evolution with the Akaike information criterion (AIC) and 548 Akaike weights. As tests for correlation of traits may be sensitive to the root state, we explored 549 the impact of applying different root states by fixing the root to represent either wingless + no 550 ocelli, winged + no ocelli, and wings + ocelli. We further tested for the correlation of ocelli and 551 either of the three wing states (apterous, micropterous, macropterous). We estimated and 552 compared the values of log-likelihoods and Akaike weights, and corrected p-values for multiple 553 comparisons using the Bonferroni correction. 554 555

Diversification analysis 556
We explored evolutionary rate dynamics with the Bayesian analysis of macroevolutionary 557 mixtures (BAMM) and the R package BAMMtools 145,146 . Priors were estimated using the 558 BAMMtools function "setBAMMpriors" and modified accordingly (expectedNumberOfShifts = 559 1.0; lambdaInitPrior = 4.80107993368211; lambdaShiftPrior = 0.00644760145800591; 560 muInitPrior = 4.80107993368211). We ran four chains of 10 million generations being sampled 561 every 1000 generations using the "speciationextinction" model and seting 562 "globalSamplingFraction" to 0.15 (~500 taxa of ~3400 described species). Convergence was 563 assessed and subsequent data analysis performed in BAMMtools. 564 Whether the evolution of wings had an impact on diversification rates was tested using the 565 Hidden State-dependent Speciation and Extinction model (HiSSE) in the R package hisse 147 . We 566 used the binary dataset of male wings and removed the outgroup (Embioptera). Following Song 567 et al. 148 , we then fitted 24 different models 147 including original BiSSE models 149 , character-568 independent diversification models (CID-2 and CID-4) and HiSSE models, which assume a 569 hidden state. All models were compared using log-likelihood, AIC and Akaike weights.     The selection of fossils to use as calibration age priors was identical for all three BI trees 1154