Butterfly wings exhibit a diversity of patterns which can vary between forewings and hindwings and spatially across the same wing. Regulation of morphological variation involves changes in how genes are expressed across different spatial scales which is driven by chromatin dynamics during development. How patterns of chromatin dynamics correspond to morphological variation remains unclear. Here we compared the chromatin landscape between forewings and hindwings and also across the proximal and distal regions of the hindwings in two butterfly species, Bicyclus anynana and Danaus plexippus. We found that the chromatin profile varied significantly between the different wing regions, however, there was no clear correspondence between the chromatin profile and the wing patterns. In some cases, wing regions with different phenotypes shared the same chromatin profile whereas those with a similar phenotype had a different profile. We also found that in the forewing, open chromatin regions were AT rich whereas those in the hindwing were GC rich. GC content also varied between the proximal and hindwing regions. These differences in GC content were also reflected in the transcription factor binding motifs that were differentially enriched between the wings and wing regions. Our results suggest that the chromatin landscape varies between different wing tissues and even spatially within the same tissue with no clear correspondence to phenotype. These findings may be explained by differences in how Hox genes cooperate with other transcription factors that show preferences for specific GC content and function either as activators or repressors in different wings or wing regions.

Wings from Bicyclus anynana and Danaus plexippus were dissected at ~22-26 hours post-pupation. For control input libraries (non-enriched), two whole forewings and two whole hindwings were pooled for each species. FAIRE-enriched libraries for D. plexippus included 3 libraries, one prepared from whole forewings, and two from partial hindwings (proximal and distal regions) (Fig.1). For B. anynana, three FAIRE-enriched libraries were prepared including a forewing distal library, and two from partial hindwings (proximal and distal regions). All FAIRE-enriched libraries were prepared from 7-8 pooled wing tissues. 

Libraries were prepared by Genotypic Technology (India), as paired-end reads (75*2) and sequenced using Illumina NextSeq. Quality checking of the raw reads was performed using FASTQC v0.11.3 and reads which had a phred score>30 were retained for downstream analyses. The reads were aligned to reference scaffolds or genomes for each species (BACs for B. anynana and the whole genome for D. plexippus) with BWA (0.7.13) using the following parameters –k INT, -w INT, -A INT, -B INT, -O INT, -E INT, -L INT, -U INT. The SAM files were converted to BAM files using SAMtools-0.1.7a and the resulting BAM files were converted to sorted BAM followed by removal of PCR duplicates. The final BAM files were then converted to BEDgraph files using BEDtools-2.14.3. Peaks were called with the MACS2 software using the aligned enriched and input (control) files with the qvalue (minimum FDR) cutoff to call significant peaks (Excel files S1-6). Fold-enrichment and log likelihood scores were calculated using the command bdgcmp script on the enriched and input BEDgraph files. The bdgcmp command also removes noise from the enriched sample relative to the control. The BEDgraph files were converted BigWig files using bdg2bw for visualization in IGV.

For this experiment we only collected one replicate per library.