Genome-Wide Gene Birth-Death Dynamics Are Associated with Diet Breadth Variation in Lepidoptera

Abstract

Comparative analyses of gene birth-death dynamics have the potential to reveal gene families that played an important role in the evolution of morphological, behavioral, or physiological variation. Here, we used whole genomes of 30 species of butterflies and moths to identify gene birth-death dynamics among the Lepidoptera that are associated with specialist or generalist feeding strategies. Our work advances this field using a uniform set of annotated proteins for all genomes, investigating associations while correcting for phylogeny, and assessing all gene families rather than a priori subsets. We discovered that the sizes of several important gene families (e.g. those associated with pesticide resistance, xenobiotic detoxification, and/or protein digestion) are significantly correlated with diet breadth. We also found 22 gene families showing significant shifts in gene birth-death dynamics at the butterfly (Papilionoidea) crown node, the most notable of which was a family of pheromone receptors that underwent a contraction potentially linked with a shift to visual-based mate recognition. Our findings highlight the importance of uniform annotations, phylogenetic corrections, and unbiased gene family analyses in generating a list of candidate genes that warrant further exploration.

Keywords: Lepidoptera; butterflies; coevolution; comparative genomics; diet breadth; gene birth–death dynamics; insect–host plant interactions; specialization.

Fig. 1.

A time-calibrated phylogeny of all species (N = 30; see supplementary table S1, Supplementary Material online for full names) used in our CAFE analyses (CAFE: compuational analysis of gene family evolution), with per-species counts of the number of host plant orders identified as being commonly used by larvae. For lists of these host plant orders, see supplementary table S8, Supplementary Material online. Non-original species silhouettes were sourced from phylopic (www.phylopic.org).

Fig. 2.

A comparison of the effects of filtering on predicted protein counts for both native and de novo (BRAKER2) annotations. The 6 letter species abbreviations here use the first 3 letters of species’ binomial nomenclature; see supplementary table S1, Supplementary Material online if necessary. A) Decrease in total protein count per species following filtering of native annotation files. Note that only 20 of 30 species included in this study had locatable native annotation files (e.g. a general feature format (GFF) file); as such, not all species have filtered native annotations. B) Decrease in total protein count per species following filtering of de novo annotation files. C) Comparison of total protein counts between filtered de novo and filtered native annotations for each species. See supplementary table S1, Supplementary Material online for full species names.

Fig. 3.

(Left) Estimated counts of genes in the 7tm odorant receptor (putative pheromone receptor) family identified as rapidly contracting at the crown of Papilionoidea. Note that all non-butterfly moths have several odorant receptor genes in this family, but most butterflies have one copy. (Right) Estimated counts of genes in the odorant/pheromone-binding protein family that was identified as rapidly contracting at the crown of Papilionoidea. For both trees, counts were taken from a root-filtered, de novo generated input where gene duplication rate was not assumed to equal gene death rate and where gene duplication and death rates were calculated independently for butterflies versus all other moths.

Fig. 4.

Effects of diet breadth on gene family size, while accounting for phylogenetic history. Diet breadth for each of n = 30 Lepidoptera species was assessed by the number of orders larvae have been observed feeding on (see supplementary table S8, Supplementary Material online). A) The effect of diet breadth is depicted as credible intervals (50%, 80%, and 95%) of posterior distributions of slopes from a multilevel model (see Methods). Each horizontal line represents a gene family, and families are ordered by their posterior averages. Results are for all gene families analyzed in the root-filtered, λ ≠ μ, BRAKER2 annotation set that had an average gene family size per species ≥ 2. The vast majority of the 1,167 gene families analyzed did not have their size significantly affected by diet breadth, as seen by their overlap with 0 (solid line). Gene families that were significantly larger in species with more host orders (generalists) or significantly larger in species with fewer host orders (specialists) are shown in the insets B) and C), respectively. MADF: myb/SANT-like domain in Adf-1. Regression lines represent fits from the multilevel model, with the shaded region representing the 95% credible interval. As many species perfectly overlap, we scaled the point size to the number of species (‘Count’) at that coordinate.