Identifying the genetic mechanisms of adaptation requires the elucidation of links between the evolution of DNA sequence, phenotype, and fitness1. Convergent evolution can be used as a guide to identify candidate mutations that underlie adaptive traits2,3,4, and new genome editing technology is facilitating functional validation of these mutations in whole organisms1,5. We combined these approaches to study a classic case of convergence in insects from six orders, including the monarch butterfly (Danaus plexippus), that have independently evolved to colonize plants that produce cardiac glycoside toxins6,7,8,9,10,11. Many of these insects evolved parallel amino acid substitutions in the α-subunit (ATPα) of the sodium pump (Na+/K+-ATPase)7,8,9,10,11, the physiological target of cardiac glycosides12. Here we describe mutational paths involving three repeatedly changing amino acid sites (111, 119 and 122) in ATPα that are associated with cardiac glycoside specialization13,14. We then performed CRISPR–Cas9 base editing on the native Atpα gene in Drosophila melanogaster flies and retraced the mutational path taken across the monarch lineage11,15. We show in vivo, in vitro and in silico that the path conferred resistance and target-site insensitivity to cardiac glycosides16, culminating in triple mutant ‘monarch flies’ that were as insensitive to cardiac glycosides as monarch butterflies. ‘Monarch flies’ retained small amounts of cardiac glycosides through metamorphosis, a trait that has been optimized in monarch butterflies to deter predators17,18,19. The order in which the substitutions evolved was explained by amelioration of antagonistic pleiotropy through epistasis13,14,20,21,22. Our study illuminates how the monarch butterfly evolved resistance to a class of plant toxins, eventually becoming unpalatable, and changing the nature of species interactions within ecological communities2,6,7,8,9,10,11,15,17,18,19.
This is a preview of subscription content
Subscription info for Chinese customers
We have a dedicated website for our Chinese customers. Please go to naturechina.com to subscribe to this journal.
Get time limited or full article access on ReadCube.
All prices are NET prices.
The data supporting the findings of this study are available within the paper and its Supplementary Information.
The code to compute RMO (reproducibility of mutational order) scores in the presence of unordered mutations and correlated pathways can be accessed in Github (https://github.com/gaguerra/ModifiedRMO). The set of R scripts implements the RMO score, first proposed by Toprak and co-workers14, with the new additions of accounting for non-independent mutational pathways (in the presence of shared ancestry) and partially unresolved mutational pathways.
Barrett, R. D. & Hoekstra, H. E. Molecular spandrels: tests of adaptation at the genetic level. Nat. Rev. Genet. 12, 767–780 (2011).
Agrawal, A. A. Toward a predictive framework for convergent evolution: integrating natural history, genetic mechanisms, and consequences for the diversity of life. Am. Nat. 190, S1–S12 (2017).
Stern, D. L. The genetic causes of convergent evolution. Nat. Rev. Genet. 14, 751–764 (2013).
Storz, J. F. Causes of molecular convergence and parallelism in protein evolution. Nat. Rev. Genet. 17, 239–250 (2016).
Siddiq, M. A., Loehlin, D. W., Montooth, K. L. & Thornton, J. W. Experimental test and refutation of a classic case of molecular adaptation in Drosophila melanogaster. Nat. Ecol. Evol. 1, 0025 (2017).
Agrawal, A. A., Petschenka, G., Bingham, R. A., Weber, M. G. & Rasmann, S. Toxic cardenolides: chemical ecology and coevolution of specialized plant–herbivore interactions. New Phytol. 194, 28–45 (2012).
Petschenka, G., Wagschal, V., von Tschirnhaus, M., Donath, A. & Dobler, S. Convergently evolved toxic secondary metabolites in plants drive the parallel molecular evolution of insect resistance. Am. Nat. 190, S29–S43 (2017).
Holzinger, F., Frick, C. & Wink, M. Molecular basis for the insensitivity of the Monarch (Danaus plexippus) to cardiac glycosides. FEBS Lett. 314, 477–480 (1992).
Dobler, S., Dalla, S., Wagschal, V. & Agrawal, A. A. Community-wide convergent evolution in insect adaptation to toxic cardenolides by substitutions in the Na,K-ATPase. Proc. Natl Acad. Sci. USA 109, 13040–13045 (2012).
Zhen, Y., Aardema, M. L., Medina, E. M., Schumer, M. & Andolfatto, P. Parallel molecular evolution in an herbivore community. Science 337, 1634–1637 (2012).
Petschenka, G. et al. Stepwise evolution of resistance to toxic cardenolides via genetic substitutions in the Na+/K+ -ATPase of milkweed butterflies (lepidoptera: Danaini). Evolution 67, 2753–2761 (2013).
Horisberger, J. D. Recent insights into the structure and mechanism of the sodium pump. Physiology (Bethesda) 19, 377–387 (2004).
Weinreich, D. M., Delaney, N. F., Depristo, M. A. & Hartl, D. L. Darwinian evolution can follow only very few mutational paths to fitter proteins. Science 312, 111–114 (2006).
Toprak, E. et al. Evolutionary paths to antibiotic resistance under dynamically sustained drug selection. Nat. Genet. 44, 101–105 (2011).
Whiteman, N. K. & Mooney, K. A. Evolutionary biology: Insects converge on resistance. Nature 489, 376–377 (2012).
Berenbaum, M. R. in Molecular Aspects of Insect–Plant Associations (eds Brattsten, L. B. & Ahmad, S.) 257–272 (Springer, 1986).
Brower, L. P., Ryerson, W. N., Coppinger, L. L. & Glazier, S. C. Ecological chemistry and the palatability spectrum. Science 161, 1349–1350 (1968).
Petschenka, G. & Agrawal, A. A. Milkweed butterfly resistance to plant toxins is linked to sequestration, not coping with a toxic diet. Proc. R. Soc. B 282, 20151865 (2015).
Zhan, S. et al. The genetics of monarch butterfly migration and warning colouration. Nature 514, 317–321 (2014).
Tarvin, R. D. et al. Interacting amino acid replacements allow poison frogs to evolve epibatidine resistance. Science 357, 1261–1266 (2017).
Hague, M. T. J. et al. Large-effect mutations generate trade-off between predatory and locomotor ability during arms race coevolution with deadly prey. Evol. Lett. 2, 406–416 (2018).
Blount, Z. D., Barrick, J. E., Davidson, C. J. & Lenski, R. E. Genomic analysis of a key innovation in an experimental Escherichia coli population. Nature 489, 513–518 (2012).
Groen, S. C. et al. Multidrug transporters and organic anion transporting polypeptides protect insects against the toxic effects of cardenolides. Insect Biochem. Mol. Biol. 81, 51–61 (2017).
Dalla, S. & Dobler, S. Gene duplications circumvent trade-offs in enzyme function: Insect adaptation to toxic host plants. Evolution 70, 2767–2777 (2016).
Lin, S., Staahl, B. T., Alla, R. K. & Doudna, J. A. Enhanced homology-directed human genome engineering by controlled timing of CRISPR/Cas9 delivery. eLife 3, e04766 (2014).
Port, F., Chen, H.-M., Lee, T. & Bullock, S. L. Optimized CRISPR/Cas tools for efficient germline and somatic genome engineering in Drosophila. Proc. Natl Acad. Sci. USA 111, E2967–E2976 (2014).
Pegueroles, C. et al. Inversions and adaptation to the plant toxin ouabain shape DNA sequence variation within and between chromosomal inversions of Drosophila subobscura. Sci. Rep. 6, 23754 (2016).
Shorrocks, B. in The Genetics and Biology of Drosophila (eds Ashburner, M. et al.) 385–428 (Academic, 1982).
Ashmore, L. J. et al. Novel mutations affecting the Na, K ATPase alpha model complex neurological diseases and implicate the sodium pump in increased longevity. Hum. Genet. 126, 431–447 (2009).
Schubiger, M., Feng, Y., Fambrough, D. M. & Palka, J. A mutation of the Drosophila sodium pump α-subunit gene results in bang-sensitive paralysis. Neuron 12, 373–381 (1994).
Trifinopoulos, J., Nguyen, L.-T., von Haeseler, A. & Minh, B. Q. W-IQ-TREE: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Res. 44, W232–W235 (2016).
Petschenka, G., Pick, C., Wagschal, V. & Dobler, S. Functional evidence for physiological mechanisms to circumvent neurotoxicity of cardenolides in an adapted and a non-adapted hawk-moth species. Proc. R. Soc. B 280, 20123089 (2013).
Paradis, E., Claude, J. & Strimmer, K. APE: Analyses of phylogenetics and evolution in R language. Bioinformatics 20, 289–290 (2004).
Pupko, T., Pe’er, I., Shamir, R. & Graur, D. A fast algorithm for joint reconstruction of ancestral amino acid sequences. Mol. Biol. Evol. 17, 890–896 (2000).
Pond, S. L. K., Frost, S. D. W. & Muse, S. V. HyPhy: hypothesis testing using phylogenies. Bioinformatics 21, 676–679 (2005).
Letunic, I. & Bork, P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 44, W242–W245 (2016).
Levy Karin, E., Ashkenazy, H., Wicke, S., Pupko, T. & Mayrose, I. TraitRateProp: a web server for the detection of trait-dependent evolutionary rate shifts in sequence sites. Nucleic Acids Res. 45, W260–W264 (2017).
Poon, A. F., Lewis, F. I., Pond, S. L. & Frost, S. D. An evolutionary-network model reveals stratified interactions in the V3 loop of the HIV-1 envelope. PLOS Comput. Biol. 3, e231 (2007).
Gratz, S. J. et al. Highly specific and efficient CRISPR/Cas9-catalyzed homology-directed repair in Drosophila. Genetics 196, 961–971 (2014).
Ponton, F., Chapuis, M. P., Pernice, M., Sword, G. A. & Simpson, S. J. Evaluation of potential reference genes for reverse transcription-qPCR studies of physiological responses in Drosophila melanogaster. J. Insect Physiol. 57, 840–850 (2011).
Taussky, H. H. & Shorr, E. A microcolorimetric method for the determination of inorganic phosphorus. J. Biol. Chem. 202, 675–685 (1953).
Petschenka, G. & Dobler, S. Target-site sensitivity in a specialized herbivore towards major toxic compounds of its host plant: the Na+K+-ATPase of the oleander hawk moth (Daphnis nerii) is highly susceptible to cardenolides. Chemoecology 19, 235–239 (2009).
Beikirch, H. Toxicity of ouabain on Drosophila melanogaster. Experientia 33, 494–495 (1977).
Ja, W. W. et al. Prandiology of Drosophila and the CAFE assay. Proc. Natl Acad. Sci. USA 104, 8253–8256 (2007).
Glass, H. W. Jr & Pan, M. L. Laboratory rearing of monarch butterflies (Lepidoptera: Danaiidae), using an artificial diet. Ann. Entomol. Soc. Am. 76, 475–476 (1983).
Malcolm, S. B. & Zalucki, M. P. Milkweed latex and cardenolide induction may resolve the lethal plant defence paradox. Entomol. Exp. Appl. 80, 193–196 (1996).
Šali, A., Potterton, L., Yuan, F., van Vlijmen, H. & Karplus, M. Evaluation of comparative protein modeling by MODELLER. Proteins 23, 318–326 (1995).
Laursen, M., Yatime, L., Nissen, P. & Fedosova, N. U. Crystal structure of the high-affinity Na+K+-ATPase-ouabain complex with Mg2+ bound in the cation binding site. Proc. Natl Acad. Sci. USA 110, 10958–10963 (2013).
Gregersen, J. L., Mattle, D., Fedosova, N. U., Nissen, P. & Reinhard, L. Isolation, crystallization and crystal structure determination of bovine kidney Na+,K+-ATPase. Acta Crystallogr. F 72, 282–287 (2016).
Janson, G., Zhang, C., Prado, M. G. & Paiardini, A. PyMod 2.0: improvements in protein sequence-structure analysis and homology modeling within PyMOL. Bioinformatics 33, 444–446 (2017).
The PyMOL molecular graphics system v.1.8 (Schrödinger LLC, 2015).
Morris, G. M. et al. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 30, 2785–2791 (2009).
Sanner, M. F. Python: a programming language for software integration and development. J. Mol. Graph. Model. 17, 57–61 (1999).
Ganetzky, B. & Wu, C. F. Indirect suppression involving behavioral mutants with altered nerve excitability in Drosophila melanogaster. Genetics 100, 597–614 (1982).
We thank V. Wagschal, who helped with the construction and testing of in vitro cell lines; E. Toprak, who provided MATLAB code; T. O’Connor, who aided in the sequestration analyses; the Essig Museum of Entomology for photographs of milkweed butterfly specimens; M. Fa and K. O’Connor-Giles for advice on the development of the fly lines; and E. LaPlante for assistance with feeding assays. D. Bachtrog, K. Mooney, P. Moorjani, M. Nachman, R. Nielsen, P. Sudmant, R. Tarvin and B. Walsh provided feedback that improved the manuscript. Access to the HPC resources of Aix-Marseille Université was supported by the Equip@Meso (ANR-10-EQPX-29-01) project of the Investissements d’Avenir supervised by the Agence Nationale de la Recherche. This project was supported by grants from the Gordon and Betty Moore Foundation (Life Sciences Research Foundation Postdoctoral Fellowship Grant GBMF2550.06 to S.C.G.), the German Research Foundation (DFG, grant Do527/5-1 to S.D.), the Agence National de la Recherche (grant no. BioHSFS ANR-15-CE11-0007 to F.S. and F.R.), the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 772257 to F.S. and F.R.), the National Geographic Society (grant 9097-12 to N.K.W.), the National Science Foundation (grant DEB-1256758 to N.K.W. and IOS-1907491 to A.A.A.), the John Templeton Foundation (grant ID 41855 to A.A.A., S.D., and N.K.W.), and the National Institute of General Medical Sciences of the National Institutes of Health (award no. R35GM119816 to N.K.W.).
The authors declare no competing interests.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Peer review information Nature thanks Joseph W. Thornton and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Extended data figures and tables
Extended Data Fig. 1 Substitutions at ATPα amino acid residues 111, 119 and 122 are directly or indirectly associated with insect specialization on plants that produce cardiac glycosides.
a, The number of occurrences of each substitution across the 21 lineages in which specialization evolved independently. b, TraitRateProp analysis of the H1–H2 loop of ATPα across insects shows amino acid residues that are strongly associated with feeding on cardiac glycoside-producing plants and toxin sequestration. Bayes factor values in the top histogram indicate per-site associations between feeding and sequence rate evolution, values in the bottom histogram indicate per-site associations between sequestration and sequence rate evolution. Values over 10 were considered different (asterisks). For information on the species included in the analysis please see Supplementary Text. Colours in the multi-sequence alignment represent individual amino acids. c, BGM shows the correlated evolution of amino acid sites within the H1–H2 loop of ATPα. The table shows the marginal posterior probabilities (PP) between amino acid interactions, where the PP exceeds a default cut-off of 0.5. The residue interactions are depicted graphically, with amino acid sites represented by the nodes and the PP associated with a given epistatic or co-evolutionary interaction indicated by the values at the arrows. Nodes circled in orange indicate amino acid sites that are the focus of experiments in this study. Sites 111 and 122 are very strongly associated with feeding and sequestering, and site 119 co-evolves with site 111.
Extended Data Fig. 2 A two-step genome editing approach using CRISPR–Cas9 and HDR to generate knock-in Atpα lines of Drosophila melanogaster.
In the first step, the region encoding the H1–H2 extracellular domain was replaced with a 3×P3::GFP marker through CRISPR–Cas9-mediated HDR. This generated a common parent line with the deletion allele AtpαDeletion (GFP+). In the second step, the 3×P3::GFP marker was replaced with each of the synonymous and non-synonymous point mutation alleles through an additional round of CRISPR–Cas9-mediated HDR to generate the knock-in lines. The crossing schemes to establish the deletion line and the knock-in lines following the first and second rounds of HDR, respectively, are also shown. See also Methods and Supplementary Tables 4–7 for further details on the genome engineering strategy and crosses behind the establishment of the knock-in lines.
Extended Data Fig. 3 Point mutations have some effect on adult emergence, but do not lead to major changes in baseline Atpα expression or Na+/K+-ATPase activity.
a, Percentages of emerging adults of the knock-in and control lines on standard Drosophila medium (n = 7–8 vials, 100 eggs per vial, mean ± s.e.m.). Survival of the knock-in lines and control lines QAN (engineered control) and QAN* (w1118 wild type) was compared using one-way ANOVA (P < 0.001). Survival differed between QAN* and some of the knock-in lines, but not between the engineered control line QAN and any of the knock-in lines except VSN (post hoc Tukey’s tests (letters)). b, Atpα expression was not different among the engineered Drosophila knock-in lines or w1118 wild-type flies (QAN*). Atpα transcript level differences were assayed by qPCR. Expression was assayed in three biological replicates (symbols represent the mean ± s.e.m.), with two technical replicates per biological replicate (averaged for each biological replicate), of five- to six-day-old females as fold change standardized against rpl32 expression in QAN* flies. The expression fold change between genotypes was compared using one-way ANOVA (P = 0.3197). c, None of the sequential Atpα genotypes found along the monarch lineage affected base-line levels of pump activity in a sodium pump enzymatic assay using extracts of fly heads (one-way ANOVA, P = 0.1377; symbols represent the mean ± s.e.m. of 3–7 biological replicates). Further information on experimental design and statistical test results is in the Source Data. ns, not significant.
Extended Data Fig. 4 Drosophila flies with edited genomes show increased larval–adult survival when fed cardiac glycosides.
These panels accompany Fig. 2c. Larval–adult survival when reared on diets with a range of ouabain concentrations was different between monarch lineage knock-in lines relative to control lines (QAN, engineered control; QAN*, w1118 wild type). Symbols represent the mean ± s.e.m. of 3–6 biological replicates (50 larvae per replicate). Curves were fit through a univariate logistic regression (effect of ouabain concentration on survival), and the difference in survivorship trajectories between each pair of fly lines (genotypes) was evaluated by performing an LRT to assess the significance of the inclusion of an interaction term between genotype and ouabain concentration in the logistic regression for a pair of lines (**P < 0.01, ***P < 0.001). Further information on experiment design and statistical test results is in the Source Data.
Extended Data Fig. 5 Drosophila flies with edited genomes show increased adult survival when fed cardiac glycosides.
These panels accompany Fig. 2d. Adult survival when reared on diets with a range of ouabain concentrations was different between monarch lineage knock-in lines and control lines (QAN, engineered control; QAN*, w1118 wild type). Symbols represent the mean ± s.e.m. of 3–6 biological replicates. Curves were fit through a univariate logistic regression (the effect of ouabain concentration on survival), and a difference in survivorship trajectories between each pair of fly lines (genotypes) was evaluated by performing an LRT to assess the significance of the inclusion of an interaction term between genotype and ouabain concentration in the logistic regression for a pair of lines (*P < 0.05, **P < 0.01, ***P < 0.001). Further information on experimental design and statistical test results is in the Source Data.
Extended Data Fig. 6 The mutational path of ATPα in the monarch butterfly lineage increases dietary tolerance to ouabain in vivo in engineered Drosophila without affecting feeding rate.
a, Estimation of mean lifespan (days) in adult females (four to seven days old at the start of the experiment) of the knock-in and control lines in CAFE assays across a range of ouabain concentrations. Each data point represents the mean ± s.e.m. of five biological replicates. Both ouabain concentration and genotype affect the survival time (two-way ANOVA (P < 0.0001) with post hoc Tukey’s tests (letters indicate pairwise differences between genotypes)). b, Estimation of LD50 (μl per fly; solid lines) and feeding rates (μl per fly per day; dashed lines) in the same individuals as in a. Each data point represents the mean ± s.e.m. of five biological replicates. Both ouabain concentration and genotype affect LD50 (two-way ANOVA (P < 0.0001) with post hoc Tukey’s tests (letters indicate pairwise differences between genotypes)). Further information on experimental design and statistical test results is in the Source Data.
Extended Data Fig. 7 Survival of knock-in lines on fly diets supplemented with dried, pulverized leaves of two milkweed species that host monarch butterflies in nature (A. curassavica and A. fascicularis).
a, Photograph of A. curassavica plant used in this study. b, c, Percentages of pupariating larvae and emerging adults of the knock-in and control (wild-type w1118: QAN*) lines on fly diet with and without A. curassavica leaf material (n = 3–4, mean ± s.e.m.). b′, c′, Differences in pupariation and emergence percentages on a fly diet with milkweed relative to percentages on a control diet (n = 3–4, mean ± s.e.m.). Mean differences between percentages in b′ and c′ were tested with one-way ANOVA (P < 0.01) followed by post hoc Tukey’s tests (letters). These panels accompany Fig. 2e. d, Photograph of A. fascicularis plant used in this study. e, f, Percentages of pupariating larvae and emerging adults of the knock-in and control lines on fly diet with and without A. fascicularis leaf material (n = 4, mean ± s.e.m.). e′, f′, Differences in pupariation and emergence percentages on a fly diet with milkweed relative to percentages on a control diet (n = 4, mean ± s.e.m.). Mean differences between percentages in e′ and f′ were tested with one-way ANOVA (P < 0.001) followed by post hoc Tukey’s tests (letters). Experiments were performed once, and adding leaf material of either of the two milkweed species to the fly diet had largely consistent effects on survival of the monarch lineage knock-in and control fly lines.
Extended Data Fig. 8 The Atpα genotypes found along the monarch lineage sequentially increase TSI to ouabain without affecting baseline levels of sodium pump activity.
a, In vitro ouabain sensitivity of engineered Drosophila Na+/K+-ATPases transiently expressed in Sf9 cell lines. Each data point represents the mean ± s.e.m. of three biological replicates. log10[IC50] for each type of Na+/K+-ATPase was estimated after four-parameter logistic curve fitting, and statistical differences between log10[IC50] values were tested with one-way ANOVA (P < 0.0001) followed by post hoc Tukey’s tests (letters). b, None of the sequential Atpα genotypes found along the monarch lineage affected baseline levels of pump activity in the enzymatic assay with extracts of transiently transfected Sf9 cells (one-way ANOVA, P = 0.3197). Each data point represents the mean ± s.e.m. of three biological replicates.
Extended Data Fig. 9 Molecular docking simulations show stepwise reductions in ouabain binding to Na+/K+-ATPases with monarch lineage substitutions in ATPα.
The ouabain binding pocket structure obtained from molecular docking simulations for each Na+/K+-ATPase with mutated ATPα. The mutated residues are shown in sticks and are labelled. The H1–H2 loop of ATPα is shown in blue. The extracellular region of the α-subunit is removed for simplicity. For the wild-type (QAN) ATPase, ouabain is shown in its co-crystal structure coordinates (white, transparent) together with its best-docked position. For all other ATPases only the best-docked positions (closest to the co-crystal structure) are shown together with ouabain’s docking position for the wild-type ATPase (dark grey). The triple-mutated VSH ATPase has two distinct docking scores: one is similar to the docking energy for the wild-type ATPase and the other has the lowest binding energy compared to all other mutated ATPases. The potential existence of both states might be related to a trend of reduced bang sensitivity for this genotype compared to some of the single-mutant genotypes. A119 is not directly part of the ouabain binding pocket, and therefore, A119S alone does not change ouabain binding. Although the consequences of A119S are relatively subtle, the mutation may disrupt the local hydrogen bonding network and cause structural or dynamic changes in the loop or in its vicinity.
About this article
Cite this article
Karageorgi, M., Groen, S.C., Sumbul, F. et al. Genome editing retraces the evolution of toxin resistance in the monarch butterfly. Nature 574, 409–412 (2019). https://doi.org/10.1038/s41586-019-1610-8
Nature Communications (2021)
Genome-wide macroevolutionary signatures of key innovations in butterflies colonizing new host plants
Nature Communications (2021)
Scientific Reports (2020)
Less Is More: a Mutation in the Chemical Defense Pathway of Erysimum cheiranthoides (Brassicaceae) Reduces Total Cardenolide Abundance but Increases Resistance to Insect Herbivores
Journal of Chemical Ecology (2020)