Evolutionary dynamics of genome size and content during the adaptive radiation of Heliconiini butterflies
This article has been Reviewed by the following groups
Listed in
- Evaluated articles (Arcadia Science)
Abstract
Heliconius butterflies, a speciose genus of Müllerian mimics, represent a classic example of an adaptive radiation that includes a range of derived dietary, life history, physiological and neural traits. However, key lineages within the genus, and across the broader Heliconiini tribe, lack genomic resources, limiting our understanding of how adaptive and neutral processes shaped genome evolution during their radiation. Here, we generate highly contiguous genome assemblies for nine Heliconiini, 29 additional reference-assembled genomes, and improve 10 existing assemblies. Altogether, we provide a dataset of annotated genomes for a total of 63 species, including 58 species within the Heliconiini tribe. We use this extensive dataset to generate a robust and dated heliconiine phylogeny, describe major patterns of introgression, explore the evolution of genome architecture, and the genomic basis of key innovations in this enigmatic group, including an assessment of the evolution of putative regulatory regions at the Heliconius stem. Our work illustrates how the increased resolution provided by such dense genomic sampling improves our power to generate and test gene-phenotype hypotheses, and precisely characterize how genomes evolve.
Article activity feed
-
-
-
-
Curiously, a contraction in the Hemocyanin superfamily was only observed in H. aoede, the only Heliconius species not to feed on pollen in our data, marking hexamerins as a potential mechanistic link to the divergent strategies for nitrogen storage in pollen-feeding Heliconius (Fig. 5A).
Interesting - this seems like a perfect complement for the analysis of Cocoonases later on! Would be very interesting to apply RELAX/aBSREL in the hypthesis-testing framework (i.e. foreground/background - discussed in later comment) to this gene family to test the hypothesis that this gene-family contraction in H. aoede is paired with an increased intensity of positive selection as compared to other Heliconius butterflies.
What sre the members of the Hemocyanin superfamily that contracted within H. aoede? Are they the same (or do they include) the two …
Curiously, a contraction in the Hemocyanin superfamily was only observed in H. aoede, the only Heliconius species not to feed on pollen in our data, marking hexamerins as a potential mechanistic link to the divergent strategies for nitrogen storage in pollen-feeding Heliconius (Fig. 5A).
Interesting - this seems like a perfect complement for the analysis of Cocoonases later on! Would be very interesting to apply RELAX/aBSREL in the hypthesis-testing framework (i.e. foreground/background - discussed in later comment) to this gene family to test the hypothesis that this gene-family contraction in H. aoede is paired with an increased intensity of positive selection as compared to other Heliconius butterflies.
What sre the members of the Hemocyanin superfamily that contracted within H. aoede? Are they the same (or do they include) the two Hemocyanins that were found to have previously undergone expansions?
-
Ancestral state reconstruction of genome size was assessed using The maximum likelihood method implemented in the R package phytools (68).
Under what model specifically was ancestral state reconstruction (ASR) conducted? The ML method in phytools can conduct ASR under three models (Brownian Motion, Ornstein–Uhlenbeck, and Early Burst). If not already done so, I might suggest fitting these alternative models of trait evolution prior to conducting ASR under the best-fit model.
Given the taxonomic scale, it might also be worth considering including fitting some multi-regime BM/OU models as implemented in OUwie: https://doi.org/10.1111/j.1558-5646.2012.01619.x.
-
All sequences from each of the four clades were realigned separately and a number of tests implemented in HyPhy (overall ω; SLAC; aBSREL; RELAX). In particular, the sign of diversifying positive selection (aBSREL) was detected by scanning all internal branches of the whole cocoonases phylogeny, correcting for multiple tests using a final P-value threshold of 0.05.
If I'm reading this correctly, it seems that for each analysis, all branches were tested for signatures of selection (with the exception of SLAC which assumes selection is constant across the entire tree)? Although RELAX and aBSREL can be used to test each branch, each are significantly more powerful when used in a hypothesis-testing framework, assigning a subset of branches to the "background" (control) and the remainder to the "foreground" (treatment), and estimating …
All sequences from each of the four clades were realigned separately and a number of tests implemented in HyPhy (overall ω; SLAC; aBSREL; RELAX). In particular, the sign of diversifying positive selection (aBSREL) was detected by scanning all internal branches of the whole cocoonases phylogeny, correcting for multiple tests using a final P-value threshold of 0.05.
If I'm reading this correctly, it seems that for each analysis, all branches were tested for signatures of selection (with the exception of SLAC which assumes selection is constant across the entire tree)? Although RELAX and aBSREL can be used to test each branch, each are significantly more powerful when used in a hypothesis-testing framework, assigning a subset of branches to the "background" (control) and the remainder to the "foreground" (treatment), and estimating parameters within each test-set. The per-branch tests are generally considered to be most useful in exploratory use-cases.
I think this use-case scenario is ideally suited to the hypothesis testing framework of these two models. Specifically, to test the hypothesis that selection on cocoonases intensified (or was relaxed) in H. aoede as compared to the rest of the Heliconius butterflies, you could treat H. aoede as the foreground, and the rest of Heliconius as the background for both RELAX and aBSREL. It's unclear whether cocoonases in Eueides or other species within Heliconiini should be experiencing the same/similar selection regime as in H. aoede, so I would suggest excluding those species from this specific analysis.
-
optimizing parameters (Romain Derelle, personal communication) for a more reliable list of single-copy orthologous groups (scOGs).
How might this approach to identify orthologs, which is tailored to the identification of single-copy orthologs, impact downstream inferences about patterns of gene-family evolution? Put another way, if the approach described here is more prone to recovering scOGs (or low-copy number OGs), then it seems likely that fewer multi-copy number (particularly large/high-copy number) OGs will be identified as a result. If so, the rate estimates of gene-family expansion/contraction using CAFE seem likely to be downwardly/upwardly biased respectively, as the majority of gene-families under consideration are particularly conserved with respect to copy-number.
I think this possibility would be worth …
optimizing parameters (Romain Derelle, personal communication) for a more reliable list of single-copy orthologous groups (scOGs).
How might this approach to identify orthologs, which is tailored to the identification of single-copy orthologs, impact downstream inferences about patterns of gene-family evolution? Put another way, if the approach described here is more prone to recovering scOGs (or low-copy number OGs), then it seems likely that fewer multi-copy number (particularly large/high-copy number) OGs will be identified as a result. If so, the rate estimates of gene-family expansion/contraction using CAFE seem likely to be downwardly/upwardly biased respectively, as the majority of gene-families under consideration are particularly conserved with respect to copy-number.
I think this possibility would be worth exploring/addressing. For instance, if the implementation of broccoli is relaxed so as to be less tailored to the recovery of scOGs, does the shape of the gene-family expansion/contraction distribution fundamentally change?
-
We illustrate how dense genomic sampling improves our resolution of gene-phenotype links, and our understanding of how genomes evolve.
Really exciting work, and an immense amount of work - congratulations! This was a fun read!
It seems that towards the conclusion of the paper, you outline how your data does not appear to be consistent with one of the prevailing hypothesized mechanisms for the digestion of pollen in Heliconius. I can't help but wonder if, using the data you present here, you might be able to point towards individual gene-families/orthogroups that are consistent with the evolution of this key-innovation?
-
bar plots show total gene counts partitioned according to their orthology profiles, from Nymphalids to lineage-restricted and clade-specific genes.
What does the distribution of mean per-species gene-count look like across all orthologs? Is there a small handful of very large orthogroups that are observed in all species and comprise the majority of the gene-counts shown here?
-
Hemocyanins are also among several GFs with evidence of divergent selection regimes (ω) between Eueides and Heliconius, alongside Trypsins, Protein kinases, P450s, Sugar transporters, Ion and ABC transporters (Fig. 5B).
What was the estimated ω for H. aoede at these gene-families as compared to those distributions estimated for Eueides and Heliconius? If involved in the evolution of pollen-feeding, I think that we would expect ω to be more similar between H. aoede and Eueides than between H. aoede and Heliconius?
-
(A, B) The contribution of TEs and CDS to genome size variation across Heliconiinae, respectively.
These two relationships should probably be assessed in a phylogenetic context, accounting for shared evolutionary history and non-independence among species (e.g. using phylogenetic least-squares regression - pgls). There is clear phylogenetic structuring in these two plots. For instance, Dione/Agraulis appear to have lower TE content for their genome sizes (shallower slope and intercept), whereas Erato seems to accumulate TEs more rapidly for their corresponding genome size as compared to other clades. Raw points and the original regression slope certainly provide taxonomic intuition, but to summarize the relationship across the entire clade, a regression slope as inferred by PGLS (for instance) should be reported.
-
We next explored the relationship between transposable elements (TEs) and genome size, and their effect on gene architecture. We found that larger genomes tend to have a distribution of intron length skewed towards longer introns (fig. S16a), with a positive correlation between median intron length and total TE content (fig. S22; Pearson’s ρ=0.72; R2=0.51). Long introns also accumulate significantly more TEs then expected by their size (fig. S23; Wilcoxon rank-sum test P value = 2.13×10-13), with the effect of changing gene structure more than gene density (fig. S22 and fig. S24).
It's difficult to contextualize these results without access to the supplementary materials (at the time I'm reading this), but it's unclear to me why the choice was made to subset introns into those that are short and long (as in Fig. 4A). I'm particularly …
We next explored the relationship between transposable elements (TEs) and genome size, and their effect on gene architecture. We found that larger genomes tend to have a distribution of intron length skewed towards longer introns (fig. S16a), with a positive correlation between median intron length and total TE content (fig. S22; Pearson’s ρ=0.72; R2=0.51). Long introns also accumulate significantly more TEs then expected by their size (fig. S23; Wilcoxon rank-sum test P value = 2.13×10-13), with the effect of changing gene structure more than gene density (fig. S22 and fig. S24).
It's difficult to contextualize these results without access to the supplementary materials (at the time I'm reading this), but it's unclear to me why the choice was made to subset introns into those that are short and long (as in Fig. 4A). I'm particularly wondering about the choice to divide introns as short/long at the median, which doesn't appear to sort them according to a natural break in the length-distribution. Dividing introns in this way, and then summing for each species and each 'TE size-class' runs the possible risk of exaggerating the observed difference in intron length and TE content.
Additionally, the relationship between TE content and total intron length (regardless of how these traits are defined) should be analyzed in a phylogenetic context (e.g. using phylogenetic least-squares regression - pgls).
-
From each OG in-paralogs were removed (custom python script RemoveInParalogFromTree.py available at https://francicco@bitbucket.org/ebablab/custum-scripts.git). If the procedure generated a single-copy OG (nscOGs) it was analysed by contrasting evolutionary pressures between Eueides and Heliconius species. The signature of selection (aBSREL) and relaxation (RELAX (114)) were performed as implemented in HyPhy.
Could removal of in-paralogs influence/impact the signatures of selection that you recover using aBSREL and RELAX? It might be worth assessing whether these results are robust to the inclusion of in-paralogs.
For instance, neofunctionalization of in-paralogs following duplication is likely to manifest signatures of positive selection (Fernández et al., 2020: https://doi.org/10.1093/molbev/msaa110, Wang et al., 2015: https://doi.or…
From each OG in-paralogs were removed (custom python script RemoveInParalogFromTree.py available at https://francicco@bitbucket.org/ebablab/custum-scripts.git). If the procedure generated a single-copy OG (nscOGs) it was analysed by contrasting evolutionary pressures between Eueides and Heliconius species. The signature of selection (aBSREL) and relaxation (RELAX (114)) were performed as implemented in HyPhy.
Could removal of in-paralogs influence/impact the signatures of selection that you recover using aBSREL and RELAX? It might be worth assessing whether these results are robust to the inclusion of in-paralogs.
For instance, neofunctionalization of in-paralogs following duplication is likely to manifest signatures of positive selection (Fernández et al., 2020: https://doi.org/10.1093/molbev/msaa110, Wang et al., 2015: https://doi.org/10.1007/s11103-015-0285-2). Alternatively, pseudogenization of in-paralogs is expected to lead to a relaxation of selection on the duplicated gene-copy (e.g. Emerling 2018: https://doi.org/10.1016/j.ympev.2017.09.016; Calderoni et al., 2016: https://doi.org/10.1038/hdy.2016.59).
-
-
-
We illustrate how dense genomic sampling improves our resolution of gene-phenotype links, and our understanding of how genomes evolve.
Really exciting work, and an immense amount of work - congratulations! This was a fun read!
It seems that towards the conclusion of the paper, you outline how your data does not appear to be consistent with one of the prevailing hypothesized mechanisms for the digestion of pollen in Heliconius. I can't help but wonder if, using the data you present here, you might be able to point towards individual gene-families/orthogroups that are consistent with the evolution of this key-innovation?
-
bar plots show total gene counts partitioned according to their orthology profiles, from Nymphalids to lineage-restricted and clade-specific genes.
What does the distribution of mean per-species gene-count look like across all orthologs? Is there a small handful of very large orthogroups that are observed in all species and comprise the majority of the gene-counts shown here?
-
From each OG in-paralogs were removed (custom python script RemoveInParalogFromTree.py available at https://francicco@bitbucket.org/ebablab/custum-scripts.git). If the procedure generated a single-copy OG (nscOGs) it was analysed by contrasting evolutionary pressures between Eueides and Heliconius species. The signature of selection (aBSREL) and relaxation (RELAX (114)) were performed as implemented in HyPhy.
Could removal of in-paralogs influence/impact the signatures of selection that you recover using aBSREL and RELAX? It might be worth assessing whether these results are robust to the inclusion of in-paralogs.
For instance, neofunctionalization of in-paralogs following duplication is likely to manifest signatures of positive selection (Fernández et al., 2020: https://doi.org/10.1093/molbev/msaa110, Wang et al., 2015: https://doi.or…
From each OG in-paralogs were removed (custom python script RemoveInParalogFromTree.py available at https://francicco@bitbucket.org/ebablab/custum-scripts.git). If the procedure generated a single-copy OG (nscOGs) it was analysed by contrasting evolutionary pressures between Eueides and Heliconius species. The signature of selection (aBSREL) and relaxation (RELAX (114)) were performed as implemented in HyPhy.
Could removal of in-paralogs influence/impact the signatures of selection that you recover using aBSREL and RELAX? It might be worth assessing whether these results are robust to the inclusion of in-paralogs.
For instance, neofunctionalization of in-paralogs following duplication is likely to manifest signatures of positive selection (Fernández et al., 2020: https://doi.org/10.1093/molbev/msaa110, Wang et al., 2015: https://doi.org/10.1007/s11103-015-0285-2). Alternatively, pseudogenization of in-paralogs is expected to lead to a relaxation of selection on the duplicated gene-copy (e.g. Emerling 2018: https://doi.org/10.1016/j.ympev.2017.09.016; Calderoni et al., 2016: https://doi.org/10.1038/hdy.2016.59).
-
All sequences from each of the four clades were realigned separately and a number of tests implemented in HyPhy (overall ω; SLAC; aBSREL; RELAX). In particular, the sign of diversifying positive selection (aBSREL) was detected by scanning all internal branches of the whole cocoonases phylogeny, correcting for multiple tests using a final P-value threshold of 0.05.
If I'm reading this correctly, it seems that for each analysis, all branches were tested for signatures of selection (with the exception of SLAC which assumes selection is constant across the entire tree)? Although RELAX and aBSREL can be used to test each branch, each are significantly more powerful when used in a hypothesis-testing framework, assigning a subset of branches to the "background" (control) and the remainder to the "foreground" (treatment), and estimating …
All sequences from each of the four clades were realigned separately and a number of tests implemented in HyPhy (overall ω; SLAC; aBSREL; RELAX). In particular, the sign of diversifying positive selection (aBSREL) was detected by scanning all internal branches of the whole cocoonases phylogeny, correcting for multiple tests using a final P-value threshold of 0.05.
If I'm reading this correctly, it seems that for each analysis, all branches were tested for signatures of selection (with the exception of SLAC which assumes selection is constant across the entire tree)? Although RELAX and aBSREL can be used to test each branch, each are significantly more powerful when used in a hypothesis-testing framework, assigning a subset of branches to the "background" (control) and the remainder to the "foreground" (treatment), and estimating parameters within each test-set. The per-branch tests are generally considered to be most useful in exploratory use-cases.
I think this use-case scenario is ideally suited to the hypothesis testing framework of these two models. Specifically, to test the hypothesis that selection on cocoonases intensified (or was relaxed) in H. aoede as compared to the rest of the Heliconius butterflies, you could treat H. aoede as the foreground, and the rest of Heliconius as the background for both RELAX and aBSREL. It's unclear whether cocoonases in Eueides or other species within Heliconiini should be experiencing the same/similar selection regime as in H. aoede, so I would suggest excluding those species from this specific analysis.
-
(A, B) The contribution of TEs and CDS to genome size variation across Heliconiinae, respectively.
These two relationships should probably be assessed in a phylogenetic context, accounting for shared evolutionary history and non-independence among species (e.g. using phylogenetic least-squares regression - pgls). There is clear phylogenetic structuring in these two plots. For instance, Dione/Agraulis appear to have lower TE content for their genome sizes (shallower slope and intercept), whereas Erato seems to accumulate TEs more rapidly for their corresponding genome size as compared to other clades. Raw points and the original regression slope certainly provide taxonomic intuition, but to summarize the relationship across the entire clade, a regression slope as inferred by PGLS (for instance) should be reported.
-
We next explored the relationship between transposable elements (TEs) and genome size, and their effect on gene architecture. We found that larger genomes tend to have a distribution of intron length skewed towards longer introns (fig. S16a), with a positive correlation between median intron length and total TE content (fig. S22; Pearson’s ρ=0.72; R2=0.51). Long introns also accumulate significantly more TEs then expected by their size (fig. S23; Wilcoxon rank-sum test P value = 2.13×10-13), with the effect of changing gene structure more than gene density (fig. S22 and fig. S24).
It's difficult to contextualize these results without access to the supplementary materials (at the time I'm reading this), but it's unclear to me why the choice was made to subset introns into those that are short and long (as in Fig. 4A). I'm particularly …
We next explored the relationship between transposable elements (TEs) and genome size, and their effect on gene architecture. We found that larger genomes tend to have a distribution of intron length skewed towards longer introns (fig. S16a), with a positive correlation between median intron length and total TE content (fig. S22; Pearson’s ρ=0.72; R2=0.51). Long introns also accumulate significantly more TEs then expected by their size (fig. S23; Wilcoxon rank-sum test P value = 2.13×10-13), with the effect of changing gene structure more than gene density (fig. S22 and fig. S24).
It's difficult to contextualize these results without access to the supplementary materials (at the time I'm reading this), but it's unclear to me why the choice was made to subset introns into those that are short and long (as in Fig. 4A). I'm particularly wondering about the choice to divide introns as short/long at the median, which doesn't appear to sort them according to a natural break in the length-distribution. Dividing introns in this way, and then summing for each species and each 'TE size-class' runs the possible risk of exaggerating the observed difference in intron length and TE content.
Additionally, the relationship between TE content and total intron length (regardless of how these traits are defined) should be analyzed in a phylogenetic context (e.g. using phylogenetic least-squares regression - pgls).
-
Hemocyanins are also among several GFs with evidence of divergent selection regimes (ω) between Eueides and Heliconius, alongside Trypsins, Protein kinases, P450s, Sugar transporters, Ion and ABC transporters (Fig. 5B).
What was the estimated ω for H. aoede at these gene-families as compared to those distributions estimated for Eueides and Heliconius? If involved in the evolution of pollen-feeding, I think that we would expect ω to be more similar between H. aoede and Eueides than between H. aoede and Heliconius?
-
Curiously, a contraction in the Hemocyanin superfamily was only observed in H. aoede, the only Heliconius species not to feed on pollen in our data, marking hexamerins as a potential mechanistic link to the divergent strategies for nitrogen storage in pollen-feeding Heliconius (Fig. 5A).
Interesting - this seems like a perfect complement for the analysis of Cocoonases later on! Would be very interesting to apply RELAX/aBSREL in the hypthesis-testing framework (i.e. foreground/background - discussed in later comment) to this gene family to test the hypothesis that this gene-family contraction in H. aoede is paired with an increased intensity of positive selection as compared to other Heliconius butterflies.
What sre the members of the Hemocyanin superfamily that contracted within H. aoede? Are they the same (or do they include) the two …
Curiously, a contraction in the Hemocyanin superfamily was only observed in H. aoede, the only Heliconius species not to feed on pollen in our data, marking hexamerins as a potential mechanistic link to the divergent strategies for nitrogen storage in pollen-feeding Heliconius (Fig. 5A).
Interesting - this seems like a perfect complement for the analysis of Cocoonases later on! Would be very interesting to apply RELAX/aBSREL in the hypthesis-testing framework (i.e. foreground/background - discussed in later comment) to this gene family to test the hypothesis that this gene-family contraction in H. aoede is paired with an increased intensity of positive selection as compared to other Heliconius butterflies.
What sre the members of the Hemocyanin superfamily that contracted within H. aoede? Are they the same (or do they include) the two Hemocyanins that were found to have previously undergone expansions?
-
optimizing parameters (Romain Derelle, personal communication) for a more reliable list of single-copy orthologous groups (scOGs).
How might this approach to identify orthologs, which is tailored to the identification of single-copy orthologs, impact downstream inferences about patterns of gene-family evolution? Put another way, if the approach described here is more prone to recovering scOGs (or low-copy number OGs), then it seems likely that fewer multi-copy number (particularly large/high-copy number) OGs will be identified as a result. If so, the rate estimates of gene-family expansion/contraction using CAFE seem likely to be downwardly/upwardly biased respectively, as the majority of gene-families under consideration are particularly conserved with respect to copy-number.
I think this possibility would be worth …
optimizing parameters (Romain Derelle, personal communication) for a more reliable list of single-copy orthologous groups (scOGs).
How might this approach to identify orthologs, which is tailored to the identification of single-copy orthologs, impact downstream inferences about patterns of gene-family evolution? Put another way, if the approach described here is more prone to recovering scOGs (or low-copy number OGs), then it seems likely that fewer multi-copy number (particularly large/high-copy number) OGs will be identified as a result. If so, the rate estimates of gene-family expansion/contraction using CAFE seem likely to be downwardly/upwardly biased respectively, as the majority of gene-families under consideration are particularly conserved with respect to copy-number.
I think this possibility would be worth exploring/addressing. For instance, if the implementation of broccoli is relaxed so as to be less tailored to the recovery of scOGs, does the shape of the gene-family expansion/contraction distribution fundamentally change?
-
Ancestral state reconstruction of genome size was assessed using The maximum likelihood method implemented in the R package phytools (68).
Under what model specifically was ancestral state reconstruction (ASR) conducted? The ML method in phytools can conduct ASR under three models (Brownian Motion, Ornstein–Uhlenbeck, and Early Burst). If not already done so, I might suggest fitting these alternative models of trait evolution prior to conducting ASR under the best-fit model.
Given the taxonomic scale, it might also be worth considering including fitting some multi-regime BM/OU models as implemented in OUwie: https://doi.org/10.1111/j.1558-5646.2012.01619.x.
-