Phylogenomics reveals the evolutionary timing and pattern of butterflies and moths.

Butterflies and moths (Lepidoptera) are one of the major superradiations of insects, comprising nearly 160,000 described extant species. As herbivores, pollinators, and prey, Lepidoptera play a fundamental role in almost every terrestrial ecosystem. Lepidoptera are also indicators of environmental change and serve as models for research on mimicry and genetics. They have been central to the development of coevolutionary hypotheses, such as butterflies with flowering plants and moths' evolutionary arms race with echolocating bats. However, these hypotheses have not been rigorously tested, because a robust lepidopteran phylogeny and timing of evolutionary novelties are lacking. To address these issues, we inferred a comprehensive phylogeny of Lepidoptera, using the largest dataset assembled for the order (2,098 orthologous protein-coding genes from transcriptomes of 186 species, representing nearly all superfamilies), and dated it with carefully evaluated synapomorphy-based fossils. The oldest members of the Lepidoptera crown group appeared in the Late Carboniferous (∼300 Ma) and fed on nonvascular land plants. Lepidoptera evolved the tube-like proboscis in the Middle Triassic (∼241 Ma), which allowed them to acquire nectar from flowering plants. This morphological innovation, along with other traits, likely promoted the extraordinary diversification of superfamily-level lepidopteran crown groups. The ancestor of butterflies was likely nocturnal, and our results indicate that butterflies became day-flying in the Late Cretaceous (∼98 Ma). Moth hearing organs arose multiple times before the evolutionary arms race between moths and bats, perhaps initially detecting a wide range of sound frequencies before being co-opted to specifically detect bat sonar. Our study provides an essential framework for future comparative studies on butterfly and moth evolution.

Butterflies and moths (Lepidoptera) are 1 of the 4 major insect superradiations, with close to 160,000 described species (1, 2). They thrive in nearly all terrestrial ecosystems, feature a wide spectrum of ecological adaptations, and are inextricably connected to the natural histories of many plants, predators, and parasitoids. The extraordinary diversity of Lepidoptera is presumed to be tightly linked to the rise of flowering plants (3), and the association between butterflies and angiosperms served as the foundation for Ehrlich and Raven's theory of coevolution (4). Prior studies of lepidopteran evolution postulated that the larvae of the earliest lineages were internal feeders on nonvascular land plants and their adults had mandibulate chewing mouthparts (5, 6). The remaining majority of Lepidoptera are hypothesized to have diversified largely in parallel with angiosperms, with larvae that are external plant feeders and adults that have tube-like proboscides (7). The proboscis promoted access to flower nectar and may have led to the ability of these ancestors to disperse to greater distances and colonize new host plants.

Nocturnal moths represent >75% of Lepidoptera species diversity (8), and they are thought to have proliferated in the Late Cretaceous when they evolved bat-detecting ultrasonic hearing organs (9, 10). At least 10 lepidopteran families have these ears (11), which are on different parts of the body in different families (). The timing of when ears arose has remained untested because a stable, dated phylogenetic framework for the Lepidoptera has been unavailable. Previous phylogenetic studies that examined macroevolutionary patterns of Lepidoptera (e.g., refs. 6 and 12–17) lacked sufficient genetic data and carefully evaluated fossil sampling to confidently resolve the evolutionary history of the order.

We constructed the largest dataset of transcriptomes for phylogenetic inference of the order, comprising 2,098 protein-coding genes sampled from 186 extant Lepidoptera species. Using different dating schemes and a set of critically evaluated fossils, we conducted multiple fossil-calibrated phylogenetic analyses. With the dated tree, we examined the evolution and timing of angiosperm feeding, ultrasonic hearing organs, and adult diel activity. Our phylogenetic analyses produced a robust reconstruction of the evolutionary history and timing of these key ecological adaptations in this megadiverse insect order.

Results and Discussion

We assembled de novo characterized transcriptomes and genomes of 186 Lepidoptera species, representing 34 superfamilies, and 17 outgroup species (Datasets S1–S4). Phylogenetic analyses were based on nucleotide and amino acid datasets of 2,098 protein coding gene alignments (2,249,363 nucleotides and 749,791 amino acid sites). In addition to maximum likelihood (ML) tree reconstruction, we conducted multispecies coalescence (MSC) analyses and 4-cluster likelihood mapping (FcLM) (18) to alleviate limitations that exist with traditional support metrics, and evaluate the presence of alternative, potentially conflicting signals in our datasets. ML analyses resulted in statistically well-supported trees (), which were largely corroborated by MSC results ( and Dataset S5). The FcLM analyses (Datasets S6–S8) indicated that there was little confounding signal. The placement of a few lineages remains uncertain ( and Dataset S9), but they do not significantly impact our conclusions. All 8 divergence time analyses resulted in largely congruent age estimates ( and Datasets S9–S11). The difference between the oldest and youngest median age estimates for crown Lepidoptera was ∼21 Ma, which is ∼7% of the median age presented in Fig. 1. When the same calculations are performed on 27 other key lepidopteran crown nodes, only 1 had a median age estimate difference >26 Ma (Dataset S11).

Fig. 1.

Dated evolutionary tree of butterfly and moth relationships. The tree is derived from a maximum-likelihood analysis of 749,791 amino acid sites. Branch lengths and node ages are computed in a time-calibrated analysis of 198,050 amino acid sites and 16 fossil calibrations. Gray bars depict 95% credibility intervals of node ages. Asterisks indicate superfamilies that are nonmonophyletic in the tree. The color-coding of nodes indicates nonparametric bootstrap support values. Two of the 16 fossils were placed on outgroup branches; their placements are shown in . Scale bar is in millions of years. MICROPT, Micropterigoidea; AGATH, Agathiphagoidea; HETEROBATH, Heterobathmioidea. Additional information on fossil calibrations is provided in Dataset S10 and .

Phylogenetic analyses, regardless of dataset or analysis type, recovered the monophyly of Lepidoptera and many other previously hypothesized clades within the order, including Angiospermivora, Glossata, and Ditrysia, with strong branch support (Dataset S5). Our results indicate that the most recent common ancestor of crown Lepidoptera appeared in the Late Carboniferous, ∼299.5 Ma (95% credibility interval [CI], 312.4 to 276.4 Ma; Fig. 1), significantly predating the oldest known lepidopteran fossil from the Triassic-Jurassic boundary (201 Ma; Dataset S10). Among Lepidoptera, the superfamilies Micropterigoidea, Agathiphagoidea, and Heterobathmioidea have been considered “primitive moths,” because their adults have functional mandibulate mouthparts and other plesiomorphic features, including feeding on nonvascular land plants (19). Prior molecular studies suggested that Agathiphagoidea, which are seed borers on the ancient conifer family Araucariaceae (19, 20), is the sister group to either Micropterigoidea (12, 21) or Heterobathmioidea (22). Our study contradicts these findings, confidently placing Micropterigoidea, which are detritivores or bryophyte feeders (7), as the sister group to all other Lepidoptera (280.6 Ma; CI, 297.0 to 257.4 Ma), with Agathiphagoidea as the closest relative to all remaining lineages, corroborating traditional morphology-based studies (19). Our results imply that the ancestral mandibulate moth may have been feeding on bryophytes ∼300 Ma, followed by some lineages transitioning to vascular plants and then to angiosperms, largely following the evolutionary history of plants.

Angiospermivora, which feed predominantly on flowering plants as larvae (21), are composed of the mandibulate Heterobathmioidea and all proboscis-bearing Lepidoptera. The crown age of Angiospermivora is dated at 257.7 Ma (CI, 276.7 to 234.5 Ma), and our analyses unequivocally inferred Heterobathmioidea as a sister group to all other taxa within this clade. Both groups have leaf-mining larvae (19), suggesting that moths were internal plant tissue feeders in the Late Permian. Although the crown age of flowering plants remains controversial (23), there is general agreement that some angiosperm lineages, which are in part now extinct, were present >300 Ma (24), a consensus corroborated by a long stem branch that links angiosperms to gymnosperms in recent studies (24–29). Our results support the hypothesis that Lepidoptera were internal feeders before becoming predominantly external herbivores (1, 3, 7, 19). Large body size, a characteristic of some extant butterflies and moths, may be an ecological consequence of being freed from the constraints of living within plants (7), thereby allowing greater access to nectar sources and the ability to fly far distances.

The adult proboscis is the unifying feature of clade Glossata, which contains >99% of all extant Lepidoptera species (2). The proboscis is considered an important innovation that enabled lepidopteran diversification (19). Adults of Heterobathmioidea, the sister group to Glossata, use their mandibles to eat pollen of the ancient Nothofagus tree (19). Adults of Eriocranioidea, part of a 3-family clade sister to all other Glossata, use their proboscis to drink water or sap (19), a behavior that may have been the precursor to nectar feeding. The common ancestor of nectar-feeding Lepidoptera first appeared in the Middle Triassic (241.4 Ma; CI, 261.1 to 218.9 Ma), overlapping with the estimated diversification periods of speciose flowering plant crown groups (Fig. 2).

Fig. 2.

Correlations among evolution of lepidopteran traits, angiosperms, and bats. Dated lepidopteran tree from Fig. 1 (right half of circle, purple tree); recent dated angiosperm phylogeny from ref. 27 (left half of circle, green tree). The blue ring represents the credibility interval for the age of bats (30). Red nodes represent independent origins of ancestral Lepidoptera with hearing organs; 4 additional origins are not shown on this tree because they are represented by terminal lineages and could not be assigned to a node (). Diurnal butterfly lineages are shown as orange branches. The large light green circle in the middle of the tree corresponds to bar 5 in the rectangular box on the right. The 8 green vertical bars in this rectangular box represent stem-to-crown intervals for the age of angiosperms based on recent evolutionary studies (Dataset S12). Tick marks on vertical scale bars represent 10-Ma intervals. Images are of the oldest known echolocating bat fossil, Icaronycteris index (A), Lepidoptera scale (B), and angiosperm pollen (C). Fig. 2A is reprinted with permission of the Royal Ontario Museum, © ROM. Fig. 2B is reprinted from ref. 66; © The Authors of ref. 66, some rights reserved; exclusive licensee American Association for the Advancement of Science. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC). Fig. 2C is reprinted with permission from ref. 69, which is licensed under CC BY 3.0. Full accreditation and additional acknowledgments are provided in .

The clade Ditrysia, which is estimated to have a late-Jurassic crown age (154.7 Ma; CI, 172.1 to 137.5 Ma) is the most taxonomically and ecologically diverse group of Lepidoptera, with >150,000 described species. Ditrysia includes butterflies and many moth superfamilies that predominantly feed on angiosperms as larvae. Monophyly of Ditrysia, recovered in all of our molecular analyses, is supported by a unique apomorphy: females have 2 reproductive openings, 1 for mating and another for egg-laying (17, 19). Relationships among major ditrysian groups such as Apoditrysia, sensu Regier et al. (12) and Mutanen et al. (22), uniting butterflies and many diverse moth families, are statistically well supported (Dataset S5). Most Apoditrysia have additional sensilla on the proboscis (19), implying that further specialization on adult food resources may have arisen in the Early Cretaceous (118.5 Ma; CI, 132.1 to 105.6 Ma).

Butterflies (Papilionoidea) are some of the most popular arthropods, and include nearly 19,000 described species (3). One of the major challenges thus far has been to confidently place butterflies within Lepidoptera (12, 17, 21, 30). We identified Papilionoidea as the sister group to all remaining Obtectomera (including silk moths, owlet moths, grass moths, and others; Fig. 1) with strong nodal support and statistical evidence from topology tests that rule out the impact of confounding signal (). This result confidently shows that butterflies are diurnal moths with a Late Cretaceous crown age (98.3 Ma; CI, 110.3 to 86.9 Ma), slightly younger than other published estimates using butterfly fossils and host-plant maximum age calibrations (13, 31–34). Swallowtails (Papilionoidea) were the sister group to the remaining butterfly families; skippers (Hesperiidae) and neotropical nocturnal moth-like butterflies (Hedylidae) were sister groups (84.4 Ma; CI, 96.6 to 72.2 Ma). Whites and sulfurs (Pieridae; 51.7 Ma; CI, 63.9 to 40.4 Ma) formed a clade that was the sister group to the clade containing the hairstreaks, blues, metalmarks (Lycaenidae + Riodinidae; 66.6 Ma; CI, 77.7 to 55.7 Ma), and brushfoot butterflies (Nymphalidae; 69.4 Ma; CI, 80.2 to 59.0 Ma).

Butterflies were hypothesized to have become diurnal because of selective pressure from predatory bats (35, 36). Phylogenetic studies that used molecular and morphological data independently confirmed the crown age of bats to be 55–65 Ma (37–41) and the crown age of the bat clade with laryngeal echolocation to be ∼50 Ma (10, 40, 42, 43). Nocturnality in Lepidoptera dates back to at least the ancestor of Heteroneura, in the Jurassic (209.7 Ma; CI, 249.9 to 208.1 Ma) ( and Dataset S11), and the switch to diurnality occurred multiple times, but most notably in the ancestor of Papilionoidea (98.3 Ma; CI, 110.3 to 86.9 Ma). Thus, diurnal activity in butterflies likely cannot be attributed to bat predation but may be due to another factor, such as nectar availability during the day. Bees, whose crown age is estimated as 125 to 100 Ma (44), may have driven angiosperm flower color evolution (36), and butterflies followed by becoming diurnal opportunists of these nectar resources.

The extraordinary diversity of ultrasonic hearing organs in nocturnal moths is thought to have evolved in response to the diversification of the echolocating-bat crown group in the Early Paleogene (11). We identified 9 different origins of hearing organs in nocturnal moth clades (), more than previously hypothesized (11). Four of these are species-rich clades (Drepanoidea, Geometroidea, Noctuoidea, and Pyraloidea), in which ears appear to have arisen in the Late Cretaceous (median age range of crown nodes, 91.6 to 77.6 Ma; CI, 103.4 to 67 Ma; Fig. 2), millions of years before echolocating bats. Our results imply that moth hearing organs might not have arisen in response to bat predation but likely evolved in response to a different selective pressure. Many moths produce sounds for sexual communication (9), but the most likely explanation for the evolution of hearing organs in Lepidoptera, and animals more broadly, is general auditory surveillance of the environment for broadband sounds produced by animal movement (e.g., predators) (45). Behavioral and neural evidence has shown that moths respond to the low-frequency (<20 kHz) walking and wingbeat sounds of predatory birds (46, 47). Day-flying Lepidoptera and moths endemic to bat-free islands show reduced auditory sensitivity to higher frequencies yet retain reactive thresholds at lower frequencies, comparable to the thresholds of other nocturnal moth species (10). Therefore, hearing organs likely evolved first for auditory surveillance before subsequently being co-opted for bat detection.

In summary, our study reveals that the common ancestor of the butterflies and moths we observe today was likely a small, Late Carboniferous species with mandibulate adults and with larvae that fed internally on nonvascular land plants. Our dating analyses with different fossils and models confirm the hypothesis that the majority of Lepidoptera proliferated with the rise of angiosperms in the Cretaceous. Hearing organs in nocturnal Lepidoptera originated several times separately, before the rise of echolocating bats. Butterflies similarly became diurnal before diversification of the bat crown group, contrary to the hypothesis that butterflies became day-fliers to escape these predators. Instead, butterflies likely became diurnal to capitalize on day-blooming flowers. Our study provides a robust phylogenetic framework and reliable time estimates for future studies on the ecology and evolution of one of the most prominent insect lineages. Data and results from this study will serve as a foundation for future comparative analyses on butterflies and moths.

Materials and Methods

Phylogenomic Analyses.

Transcriptome and genome sequences from 186 species of Lepidoptera were used for our study, along with sequences of 17 outgroup species representing other holometabolous insect orders, totaling 203 species (Dataset S1). Transcriptomes from 69 of these species were newly sequenced by either the 1K Insect Transcriptome Evolution (1KITE) consortium or by the Kawahara Lab, Florida Museum of Natural History (FLMNH). RNA extractions, mRNA isolation, fragmentation, cDNA library construction, and transcriptome sequencing were performed using the protocols of Misof et al. (48) and Peters et al. (44) for the 1KITE consortium samples. Protocols of Kawahara and Breinholt (16) were used for the FLMNH samples ().

Transcriptome assembly and contaminant removal protocols largely followed Peters et al. (44) and Kawahara and Breinholt (16); additional information on these methods is available in . An ortholog set was compiled based on the OrthoDB v7 database (49), using 2 ingroups (Bombyx mori and Danaus plexippus) and 1 outgroup (Tribolium castaneum) as reference species. The assemblies were searched for transcripts of 3,429 single-copy protein coding genes using Orthograph beta4 (50). Orthograph identified a total of 3,427 of the 3,429 orthologs across all included taxa. Amino acid sequences of each ortholog were aligned in MAFFT v7.294 (51) and refined using the protocols of Misof et al. (48). Corresponding nucleotide alignments were generated using a modified version of Pal2Nal v14 (52) with amino acid alignments as blueprints. Supermatrices for the amino acid and nucleotide datasets were generated using FasConCat-G v1.0 (53), and synonymous signal was removed from the nucleotide dataset using Degen v1.4 (54). Additional information on generation of the ortholog set and the final amino acid and “degen1” nucleotide datasets is available in . PartitionFinder2 v2.1.1 (55) and RAxML v8.2.11 (56) were used to obtain the optimal partitioning scheme for the concatenated amino acid dataset; the boundaries of the merged partitions in this scheme were adjusted at the nucleotide level and used for the concatenated degen1 nucleotide dataset. The best-fitting models were reestimated using ModelFinder in IQ-TREE v1.5.5 (57, 58) for both datasets, and maximum likelihood (ML) tree inferences were subsequently conducted. A minimum of 50 ML tree searches were performed for each dataset using IQ-TREE v1.5.5 (57); the tree with the best log-likelihood was selected as the best tree for each analysis. Multiple metrics of statistical support, including nonparametric bootstrap replicates, SH-aLRT (59), and TBE support values (60), were generated to assess the reliability of the best ML tree in each analysis. Multispecies coalescent analyses were also performed on the amino acid dataset in ASTRAL-III v5.6.3 (61), using 2 sets of individual gene trees that had been inferred with IQ-TREE v1.6.10 (57): a set of the best ML trees for each gene (inferred from 25 ML tree searches per locus) and a set of the consensus trees for each gene partition, with each consensus tree generated using 1,000 ultrafast bootstrap (UFBoot2) replicates (62).

The placement of 3 selected Lepidoptera clades (Alucitidae + Pterophoridae, Bombycoidea, and Papilionoidea) was further assessed using FcLM in IQ-TREE v1.6.7 (57), with predefined groups of taxa, to check for hidden signal and bias due to nonrandomly distributed missing data and (among-lineage) compositional heterogeneity, following the protocol of Misof et al. (48). Additional information on data partitioning, model selection, topology tests, tree inference using the partitioned supermatrix and the multispecies coalescent approach, and branch support is available in .

Divergence Time Estimation.

Divergence time estimations were computed on truncated datasets to allow completion of the analyses. AliStat v1.6 (63) was used to generate a subsampled dataset containing 195 Amphiesmenoptera species (all non-Amphiesmenoptera removed) and only including sites for which at least 80% of samples had unambiguous amino acids. MCMCTree and codeml (both part of the PAML software package, v4.9g) (64) were used to estimate divergence dates on this subsampled dataset. The best ML tree inferred from the concatenated amino acid dataset was used as the input tree. The input tree was first calibrated using age estimates of 16 carefully selected fossils, following the best-practice recommendations by Parham et al. (65). Although all 16 fossils have diagnostic morphological characteristics that enable reasonably confident placement on the tree, only 3 of these fossils have true synapomorphies. To strictly follow the guidelines of Parham et al. (65), additional analyses were performed using only these 3 fossils. We applied a conservative age constraint on the root of the input tree, with a minimum age of 201 Ma, based on the stem Glossata scale fossils discovered by van Eldijk et al. (66) and a maximum age of 314.4 Ma, based on the absence of Amphiesmenoptera fossils in the Late Carboniferous. We used 2 well-established approaches to convert fossil ages into calibrations on tree nodes. For the conservative strategy, fossil calibrations were treated with uniform distributions constrained between the corresponding fossil age (the minimum bound) and a hard maximum equal to the maximum root age. For the second strategy, the truncated-Cauchy distribution (67) was used to set calibrations for internal nodes younger than 80 Ma. (Additional information on this approach is provided in .) We applied both uncorrelated rates and autocorrelated rates to estimate divergence, for a total of 8 analyses: 2 sets of fossil calibrations (16 fossils and 3 fossils) × 2 fossil calibration strategies (uniform and Cauchy priors for nodes <80 Ma) × 2 rate types (independent and autocorrelated). To compare the radiation of Lepidoptera and flowering plants, estimates of the mean age of the ancestral angiosperm were compiled from the literature; these angiosperm ages are presented in Dataset S12 and shown in Fig. 2. Since the estimated timespan between divergence of angiosperms and gymnosperms is large, and it is possible that flowering plants existed long before the crown of angiosperms, we present the interval between the mean age of the crown (node Angiosperm) and the mean age of the stem (node Angiosperm + Gymnosperm) from the abovementioned studies (Dataset S12). Additional information and justification for our approaches to divergence time estimation is available in .

Ancestral State Reconstruction.

We conducted ancestral state reconstruction on ultrasonic hearing organs and adult diel activity. All ancestral state reconstruction analyses were performed on the best tree from the ML analysis of the concatenated amino acid dataset (), using stochastic character mapping, with the “make.simmap” command in the R package Phytools v06-44 (68). Ten thousand stochastic maps were generated for each analysis. The presence of a hearing organ was treated as a binary character in one ASR analysis and as a 7-state character in a separate ASR analysis that accounts for variation in the location of the hearing organ (see matrices A and B in Dataset S13). Diel activity was treated as a 4-state character (Dataset S14). Additional information on ancestral state reconstruction and how characters were coded is provided in .

Data Availability.

Additional data are provided in Supplementary Archives on DRYAD (70).