Genome-Wide Regulatory Adaptation Shapes Population-Level Genomic Landscapes in Heliconius

Abstract

Cis-regulatory evolution is an important engine of organismal diversification. Although recent studies have looked at genomic patterns of regulatory evolution between species, we still have a poor understanding of the magnitude and nature of regulatory variation within species. Here, we examine the evolution of regulatory element activity over wing development in three Heliconius erato butterfly populations to determine how regulatory variation is associated with population structure. We show that intraspecific divergence in chromatin accessibility and regulatory activity is abundant, and that regulatory variants are spatially clustered in the genome. Regions with strong population structure are highly enriched for regulatory variants, and enrichment patterns are associated with developmental stage and gene expression. We also found that variable regulatory elements are particularly enriched in species-specific genomic regions and long interspersed nuclear elements. Our findings suggest that genome-wide selection on chromatin accessibility and regulatory activity is an important force driving patterns of genomic divergence within Heliconius species. This work also provides a resource for the study of gene regulatory evolution in H. erato and other heliconiine butterflies.

Fig. 1.

Heliconius populations, experimental design, and data overview. (A) Three populations were used in this study: Heliconius erato petiverana (Costa Rica), H. e. lativitta (Ecuador), and incipient species H. himera (Ecuador). Black arrow indicates direct introgression, and gray arrows indicate indirect introgression. Approximate divergence for each population pair derived from Van Belleghem et al. (2017) based on nucleotide sequence. (B) Pruned phylogenetic tree based on Van Belleghem et al. (2017) shows evolutionary relationships between the study populations. (C) Experimental design for this study. (D) ChIP-seq tracks for mid-pupal forewings showing normalized read depth for H3K27ac (top, oranges) and H3K4me3 (bottom, purples) across ∼900 kb of chromosome 19. p, H. e. petiverana; l, H. e. lativitta; h, H. himera. (E) Enrichment profiles for combined H3K27ac and H3K4me3 mid-pupal histone ChIP-seq (top) and ATAC-seq (bottom) tracks show a high degree of similarity between ATAC-seq and ChIP-seq normalized read depth across all loci. This is confirmed in (F), with 91% of H3K27ac and 95% of H3K4me3 ChIP-seq peaks overlapping one or more assayed ATAC-seq peaks by at least 1 bp (aggregate of all samples and stages shown). Correlation plots of Log₁₀ signal in forewings and hindwings of H. e. lativitta show a high degree of similarity between wing tissues at both mid-pupal (G) and late-pupal (H) stages.

Fig. 2.

Variation in regulatory loci between three populations of Heliconius erato. Examples of VREs displaying population-level variability in normalized read depth for (A) ATAC-seq, (B) H3K27ac ChIP-seq, and (C) H3K4me3 ChIP-seq. (D) Percent of all unique regulatory marks displaying variability between one or more populations by assay. (E) Boxplots showing fold change between population signals for all VRE comparisons. Outliers removed due to very high fold change in the top 1% of data sets representing complete gain or loss of peaks between populations. (F) Distance of VREs from the nearest annotated TSS for ATAC-seq, H3K27ac, and H3K4me3 marked loci. (G) Distance to the next closest VRE for VREs in mid-pupal H. e. lativitta wing tissue compared with the same peak sets randomly shuffled throughout the genome (dark blue boxplots). Significance calculated by two-sample Kolmogorov–Smirnov test. For all panels: p, H. e. petiverana; l, H. e. lativitta; h, H. himera.

Fig. 3.

Genome-wide evidence of selection around VREs. (A) Genome-wide F_st for all 5-kb bins (marked “G”) and for the top 5% outlier bins (marked “O”) for all three population comparisons shows much greater population structure in outlier bins. (B) Percent of active outlier bins (showing any regulatory activity) with VREs for all three population comparisons and assays shows significant overrepresentation in outlier-associated late-pupal wings relative to mid-pupal wings. Significance determined by Fisher’s exact test. (C) Fold enrichment of VREs in outlier bins relative to active outlier genome size shows greater concentration of high F_st VREs in mid-pupal wings over late-pupal wings. (D) Percent of high F_st VREs that overlap with the bottom 2.5% of Tajima’s D bins in either population shows large enrichment of loci under strong directional selection relative to a random expectation. Significance of differences between mid-pupal and late-pupal samples in both (C) and (D) calculated by Student’s t-test. (E) Selected enriched transcription factor binding motifs within TAS variants in mid-pupal wings suggest selection on structural, olfactory, and neural phenotypes (see supplementary table S1, Supplementary Material online). For all panels: p, Heliconius erato petiverana; l, H. e. lativitta; h, H. himera.

Fig. 4.

Population structure of VRE clusters. (A) VREs in each population and developmental stage, assigned by the population showing greatest signal in pairwise population comparisons for each assay. Patterns of shared VREs between populations for (B) ATAC-seq, and (C) H3K27ac and (D) H3K4me3 ChIP-seq loci show phylogenetic distance as a major determinant of shared VREs. Gray areas indicate that the union set does not exist. Percent values show percent of elements relative to the population with the fewest assigned VREs. Percent of VREs clustering between populations, defined as nonoverlapping VREs with one or more additional VREs in another population <5 kb away, in mid-pupal and late-pupal developmental stages for (E) ATAC-seq, (F) H3K27ac, and (G) H3K4me3 ChIP-seq assays. Double bars are shown for each population pair, corresponding to the order on the pair label. Clustered VREs were categorized by evidence of either epistatic coevolution (H, top) or parallel evolution between populations (H, bottom). Gray outline shows distal clustered VREs and purple marks the putative promoter, which appears to track with changes in neighboring distal activity. The percent of clustered sites showing evidence of epistatic coevolution decreases greatly with increased phylogenetic divergence for ATAC-seq (I), and both ChIP-seq marks (J, K). The opposite trend was observed for parallel evolution in all three data types, which becomes prevalent as population divergence increases (L–N). For all panels: p, Heliconius erato petiverana; l, H. e. lativitta; h, H. himera.

Fig. 5.

Genomic features associated with regulatory variation in Heliconius. (A) Fraction of variable and nonvariable (“stable”) regulatory loci overlapping one or more SNPs. “Random” represents randomly selected genomic sequences overlapping one or more SNPs. Stable loci show a significant increase in SNPs relative to variable regulatory loci (chi-square test). (B) Conservation of VRE DNA within the species (Heliconius erato demophoon) and genus (H. melpomene) shows that most VREs are found in DNA conserved within species, but often not conserved between species. (C) Enrichment of VREs in repetitive sequences. VREs are significantly enriched in repetitive sequences, especially for variable histone marks. Bar graphs show the trend of conserved repeat associated VREs between recently diverged populations. (D) Four long interspersed nuclear element types were found to be enriched over expected frequency for VREs in the reference population. For all panels: p, H. e. petiverana; l, H. e. lativitta; h, H. himera.

Fig. 6.

Regulatory variants associate with two types of gene expression changes. (A) Example of variable regulatory activity associated with differences in gene expression between Heliconius erato petiverana and H. himera on chromosome 15 (normalized read depth shown). (B) Venn diagrams showing the number of highly differentially expressed genes in mid-pupal wing tissue genes with one or more VREs within 25 kb of the TSS. Intersections show that 82–95% of highly differentially expressed genes have a proximal VRE. (C) The majority of VREs proximal to differentially expressed genes are TAS elements, followed by H3K4me3 and H3K27ac variants, and VREs show greater association with highly differentially expressed genes compared with all genes. H3K27ac and H3K4me3 are associated with distinct changes in gene expression, with (D) H3K27ac variants associated with greater expression fold change. (E) H3K4me3 VRE-associated highly differentially expressed genes show reduced expression-normalized variance (coefficient of variation) within populations. Significance calculated by two-sample Kolmogorov–Smirnov test. For all panels: p, H. e. petiverana; l, H. e. lativitta; h, H. himera.