Butterfly eyespots evolved via cooption of an ancestral gene-regulatory network that also patterns antennae, legs, and wings

Butterfly Husbandry.

B. anynana were maintained in laboratory populations and reared at 27 °C and 60% humidity inside a climate room with a 12:12-h light:dark cycle. All larvae were supplied with young corn leaves to complete their development until pupation. Following pupation, the pupae were collected and placed in a separate cage until they emerged. The butterflies were fed every other day with banana on moist cotton in Petri dishes.

Wing Library Preparation and FAIRE-Seq Analysis.

Wings were dissected from B. anynana at ∼22 to 26 h postpupation. For control input libraries (nonenriched), two whole forewings and two whole hindwings were pooled. Three FAIRE-enriched libraries were prepared in total, including a forewing distal library (the pupal wing was cut in half and the distal region was used for the library) and two hindwing libraries, using both the proximal and distal regions of the wing. All FAIRE-enriched libraries were prepared from seven to eight pooled wing tissues. Libraries were prepared by Genotypic Technology (India) as 75-bp pair-end reads and sequenced using Illumina NextSeq. Raw reads were quality checked and reads with phred scores >30 were retained for downstream analyses. Following the removal of adapters and low-quality bases, the reads were aligned to a B. anynana BAC sequence containing Dll, with Burrows-Wheeler Aligner (0.7.13) (30), using the following parameters: -k Integer (INT), -w INT, -A INT, -B INT, -O INT, -E INT, -L INT, -U INT. The resulting SAM files were converted to binary alignment map (BAM) files, using SAMtools-0.1.7a (31). The BAM files were converted to sorted BAM, followed by removal of PCR duplicates. The final BAM files were converted to BEDgraph files, using BEDtools-2.14.3 (32). Peaks were called with MACS2 software (33), using the aligned enriched and input (control) files with the q value (minimum false discovery rate [FDR]) cutoff to call significant regions. The command bdgcmp script was used on the enriched and input BEDgraph files to generate fold enrichment and log likelihood scores. This command also removed noise from the enriched sample relative to the control. The BEDgraph files were converted to BigWig files for visualization in Integrative Genomic Viewer (IGV).

Identifying CREs for CRISPR-Cas9 Experiments.

The FAIRE-seq data were visualized using IGV. All 18 candidate CREs identified around the Dll locus were blasted against the B. anynana genome in LepBase to verify whether they were unique in the genome. Most of the candidate CREs were not unique and had multiple hits throughout the genome. One of the unique regions, the CRE Dll319, was selected as a suitable target for CRISPR knockout.

Single Guide RNA Design and Production.

sgRNA target sequences for sal were selected based on their GC content (around 60%) and the number of mismatch sequences relative to other sequences in the genome (more than three sites). In addition, we selected target sequences that started with a guanidine for subsequent in vitro transcription by T7 RNA polymerase. sgRNA for the Dll319 CRE were designed using CRISPR Direct (34), corresponding to GGN20NGG. We designed four guides spanning the length of the CRE (Fig. 3B and SI Appendix, Fig. S18 and Table S2). Two guides were designed targeting the sal740 region (SI Appendix, Fig. S28 and Table S2). The sgRNA templates were created by a PCR with overlapping primers, using Q5 polymerase (New England Biolabs [NEB]). Constructs were transcribed using T7 polymerase and (10×) transcription buffer (New England Biolabs), RNase inhibitor (Ribolock), nucleoside triphosphates (NTPs) (10 mM), and 600 ng of the PCR template. The final sample volume was 40 μL. Samples were incubated for 16 h at 37 °C and then treated with 2 μL of DNase 1 at 37 °C for 15 min. Samples were purified by ethanol precipitation, and RNA size and integrity were confirmed by gel electrophoresis.

Cas9 mRNA Production.

The plasmid pT3TS-nCas9n (Addgene) was linearized with XbaI (NEB) and purified by phenol/chloroform purification and ethanol precipitation. pT3TS-nCas9n was a gift from Wenbiao Chen, Vanderbilt University, Nashville, TN (Addgene plasmid #46757; addgene.org/46757; Research Resource identifiers [RRID]:Addgene_46757). In vitro transcription of mRNA was performed using the mMESSAGE mMACHINE T3 kit (Ambion). One microgram of linearized plasmid was used as a template, and a poly(A) tail was added to the synthesized mRNA by using the Poly(A) Tailing Kit (Thermo Fisher). The A-tailed RNA was purified by lithium-chloride precipitation and then dissolved in RNase-free water and stored at −80 °C. The Cas9 transcript was used for producing sal crispants and for the analysis of regulatory interactions among Dll, Antp, and sal.

In Vitro Cleavage Assay for the Dll319 CRE.

The sgRNAs were tested using an in vitro cleavage assay. Wild-type genomic DNA was amplified using primers that were at least 200 bp from the sgRNA sites. sgRNA (200 ng/μL per guide), Cas9 protein (800 ng/μL) (stored in a buffer containing 300 mM NaCl, 0.1 mM edetate disodium salt dihydrate, 1 mM dithiothreitol, 10 mM Tris⋅HCl, 50% glycerol pH 7.4 at 25 °C), NEB buffer 3 (1 μL), and bovine serum albumin (BSA) (1 μL) were brought to a final volume of 10 μL with nuclease-free water and incubated at 37 °C. After 15 min of incubation, the purified amplicon (100 ng) was added to the sample, which was then incubated for an additional 1 to 2 h at 37 °C. The entire reaction volume was analyzed on a 1% agarose gel. Cas9 protein was purchased from NEB EnGen Cas9 nuclear localization sequence. The cleavage assay confirmed that each guide successfully cleaved the PCR amplicon.

Embryo Injections.

Wild-type laboratory populations of B. anynana adults were provided with corn plants to lay eggs. The eggs were collected within 1.5 h of oviposition and placed onto 1-mm-wide strips of double-sided tape in plastic Petri dishes (90 mm). Cas9 protein (final concentration 800 ng/μL) and sgRNA (final concentration 200 ng/μL per guide) for all four guides were prepared in a total volume of 10 µL and incubated for 15 min at 37 °C prior to injection along with 0.5 µL of food dye to improve visualization of the injected sample into the embryos. For sal crispants, Cas9 mRNA (500 μg/μL final concentration) and sgRNA (500 μg/μL final concentration) were injected along with one-tenth of the volume of food dye. For sal740 CRE, eggs were injected with the mix of Cas9 protein (final concentration 800 ng/μL) and sgRNA (final concentration 400 ng/μL per guide). The injection mixture was kept on ice after the incubation and prior to injection. Embryo injections were carried out by nitrogen-driven injections through glass capillary needles. Injected eggs were stored in closed Petri dishes containing cotton balls that were dampened daily to maintain humidity. The hatched larvae were reared in small paper cups for 1 wk and then moved to corn plants to complete their development. SI Appendix, Tables S1, S3, and S4 summarize the injection results.

In Vivo Cleavage Assay and Genotyping of sal Crispants.

Genomic DNA was extracted with a sodium dodecyl sulfate (SDS)-based method from a pool of five injected embryos that did not hatch. About 250 bp of sequence spanning the target sequence was amplified with PCRBIO Taq Mix Red (PCR Biosystems), and PCR conditions were optimized until there were no smears, primer dimers, or extra bands. Primers are listed in SI Appendix, Table S2. The PCR products were purified with the Gene JET PCR purification kit (Thermo Fisher). A total of 200 ng of PCR product was denatured and reannealed in 10× NEB2 buffer. One microliter of T7 endonuclease I (NEB) was added to the sample, while 1 µL of MQ water was added to a negative control. Immediately after the incubation for 15 min at 37 °C, all the reactions were analyzed on a 3% agarose gel. Amplicons that showed positive cleavage from the T7 endonuclease I assay were subcloned into the pGEM-Teasy vector (Promega) through Thymine and Adenine (TA) cloning. For each target, we picked eight colonies, extracted the plasmid with a traditional alkali-SDS method, and performed a polyethylene glycol (PEG) precipitation. Sequence analysis was performed with the BIGDYE terminator kit and a 3730xl DNA Analyzer (Thermo Fisher).

Screening and Genotyping Dll319 Crispants.

Newly emerged caterpillars were screened under a microscope to look for developmental defects affecting any regions where Dll is expressed, such as the thoracic legs, mouthparts, and prolegs. Any caterpillars exhibiting defects were imaged and reared individually in paper cups until the butterflies eclosed. Caterpillars that died were immediately frozen for DNA isolation and genotyping. All other surviving caterpillars with no apparent developmental abnormalities were reared in groups on corn plants and fed ad libitum every 2 d until pupation. The eclosed butterflies were frozen individually at −20 °C. Each butterfly was carefully screened under a microscope and examined for asymmetric crispant phenotypes, focusing particularly on phenotypes expected for a Dll knockout, such as appendage, eyespot, or pigmentation defects.

Colony PCR to Identify CRE Deletions.

For selected crispants, genomic DNA was extracted from the thorax (E.Z.N.A. tissue DNA kit) and used for PCR to prepare samples for genotyping. The samples were visualized on a gel to confirm the correct size band and the PCR product was purified using a Thermo Scientific PCR purification kit. The DNA was cloned into a pGEM T-Easy Vector (Promega) and the plasmid was transformed into DH5 alpha Escherichia coli. White colony selection was used for colony PCR. The bands were visualized on a 1% agarose gel to look for bands with shifts relative to the WT band. PCR products from colonies showing evidence of a deletion were submitted for Sanger sequencing PCR (Axil Scientific, Singapore), including a sample that was amplified from B. anynana wild-type genomic DNA.

Butterfly Enhancer Reporter Assay.

A 917-bp region containing the Dll319 CRE was cloned into the piggyGUE vector via Gateway technology (Thermo Fisher). piggyGUE is the EGFP version of piggyGUM, the piggyBac-based reporter construct that was previously published (24). The details of piggyGUE will be published elsewhere. The 917-bp region was amplified from B. anynana wild-type genomic DNA using a primer containing CACC at the 5′ end for directional cloning. The PCR product was cloned into the pENTR vector and further cloned into the piggyGUE vector via a ligation reaction, as described by Lai et al. (24). A total of 4 μL of the LR reaction mix was used for bacterial transformation. After sequence analysis to confirm the presence of SNPs in the Dll319 CRE, plasmid DNA was amplified, using a Midiprep kit (Qiagen). The piggyGUE Dll319 CRE plasmid was diluted to 1 µg/µL and mixed in a 1:1 ratio with a hyperactive piggyBac transposase plasmid (35). Embryos (n = 550) were collected from B. anynana butterflies reared at 27 °C and injected ∼1 h after egg laying with the plasmid solution and a small amount of food dye, using a glass injection needle and nitrogen gas pressure. Eggs were transferred in a Petri dish to a chamber and kept moist to prevent dehydration. From this batch of eggs, 40 caterpillars hatched and were reared in paper cups during the first week and then transferred to cages with corn plants to complete their development. At all stages, caterpillars were fed corn ad libitum. From this batch of caterpillars, 19 reached adulthood (10 females and 9 males). These butterflies were evenly distributed into four cages (∼5/cage) and placed with respective wild-type males and females for breeding. We were unable to observe any dsRed signal (the positive marker of transgenesis driven by the 3xP3 promoter) in the eyes of the caterpillars from the F1 or F2 generation, despite ubiquitous dsRed signal in some first instar larvae (only) of the F1 generation, which were used later for outcrossing to wild-type individuals. This ubiquitous signal was not observed again in the offspring of these larvae. We collected eggs from the F3 generation and dissected some embryos for EGFP antibody staining. Two of the four dissected embryos showed expression of EGFP driven by the Dll319 CRE in the embryonic antennae, mouthparts, thoracic legs, and pleuropodia (Fig. 3K and SI Appendix, Fig. S23). Subsequent hemolymph PCR genetic screening in individuals of the fourth generation failed to identify additional positive individuals and the line was lost. A second line was created, using similar methods as described above, where we were able to observe EGFP expression in the eyespot centers at both larval and pupal stages (n = 6) and also in embryos (n = 2).

Drosophila Enhancer Reporter Assay.

The same 917-bp sequence that contained the Dll319 CRE was directionally cloned into pENTR-D, then Gateway cloned into the piggyPhiGUGd, the Gal4-delta version of the previously reported piggyBac-based reporter construct (24). piggyPhiGUGd also has an attB site, allowing phiC31 transgenesis. For Drosophila transgenesis, the piggyPhiGUGd Dll319 CRE construct was transformed into the attP2 site (68A4) through phiC31 integrase-mediated transgenesis system with EGFP as a visible marker (BestGene Drosophila transgenic service). Established transgenic flies were crossed with G-TRACE (36) to visualize the tissues with CRE activities.

Antibody Staining of B. anynana Embryos and Wings.

Two-day-old embryos, as well as fifth instar larval and pupal wing tissues were dissected in phosphate buffered saline (PBS) buffer under the microscope. The samples were fixed in 4% formaldehyde/fix buffer (0.1 M Pipes pH 6.9, 1 mM ethylene glycol tetraacetic acid pH 6.9, 1.0% Triton X-100, 2 mM MgSO₄) for 30 min on ice. The samples were washed with 0.02% PBSTx (PBS + Triton X-100) three times every 10 min and then blocked in 5% BSA/PBSTx for 1 h. The samples were then incubated in 5% BSA/PBSTx with the primary antibody, and incubated at 4 °C overnight. As primary antibodies, we used a rabbit polyclonal anti-Dll (at 1:200, a gift from Grace Boekhoff-Falk, University of Wisconsin, Madison, WI), a mouse monoclonal anti-Antp 4C3 (at 1:200; Developmental Studies Hybridoma Bank), a rabbit anti-Sal (at 1:20,000 for wings and pupal tissues, and 1:2,000 for embryos) (37), and a rabbit anti-EGFP antibody (at 1:200; Abcam ab290) for the transgenic embryos at 24 h (n = 4), larval wing discs from the fifth instar (n = 3), and pupal wings at 24 h (n = 3), as well as WT controls (n = 2). For double staining, we added two primary antibodies to the same tube. The wings were washed with PBSTx three times every 10 min. Then, we replaced the PBSTx with 5% BSA/PBSTx to block for 1 h, followed by incubation with the secondary antibody (1:200) in 5% BSA/PBSTx at 4 °C for 2 h. The wings were washed with PBSTx three times every 10 min, followed by mounting the wings in ProLong Gold Antifade Mountant (Thermo Fisher). The images were taken under an Olympus FV3000 confocal laser scanning microscope.

Sample Collection and Library Preparation for RNA Sequencing.

In order to identify gene expression patterns specific to eyespot formation on the developing wings, we extracted RNA from 16 different tissue types: 3- to 4-h-old, 12- to 13-h-old, and 24- to 25-h-old embryos; T1 legs, prolegs, forewings, and horns from wandering caterpillars; T1 legs, antennae, forewings, forewing margins, eyes, eyespots, and two control tissues adjacent to eyespot centers from 3-h-old pupae (Fig. 1A). For wing wounding experiments, we poked one wing between 17 and 18 h after pupation in two different places in the M3 sector, using a fine tungsten needle with a diameter of 0.25 mm and 0.001 mm at the tip (FST No. 10130-10). We collected the wings 6 h later, which corresponds to 23 to 24 h after pupation, and we also collected the other nonpoked wings as controls (5). In addition, we also collected tissues from 24-h-old and 48-h-old forewings for the clustering analysis. We performed the experiments with four biological replicates for each tissue type with 10 to 25 female individuals in each replicate (both left and right tissues were used, except for the wounded pupal wings, where a single wing was used) (SI Appendix, Table S5). Total RNA was extracted in 70 µL of nuclease-free water, using the Qiagen RNA Plus Mini kit. RNA quantity and integrity were measured using a Nanodrop and an RNA bleach gel (38). RNA libraries were prepared using the Truseq stranded mRNA kit from Illumina. Forty million reads were sequenced for each replicate using NovosEq. 6000, with 150-bp paired end and an average insert size of 250 to 300 bp. Library preparation and sequencing were carried out at AIT Novogene, Singapore. In order to avoid batch effects, we randomized the sample extraction and RNA isolation, such that two replicates of the same group were never extracted at the same time.

RNA-Seq Analysis.

The raw RNA-seq data were quality controlled and filtered. Adapter sequences and reads with low quality (less than Q30) were trimmed, using bbduk scripts (ktrim = r, k = 23, mink = 11, hdist = 1, tpe, tbo, qtrim = rl, trimq = 30, minlen = 40). In order to remove any bacterial contamination in the samples, we used the bbsplit script, which is a part of the bbmap tools (39). All bacterial genomes were downloaded from the National Center for Biotechnology Information (NCBI) (last downloaded in June 2018), and the reads were mapped to the bacterial genomes, using bbmap. Only reads whose pairs also passed through the filter were further analyzed. To remove any ribosomal RNA sequences from the RNA-seq data, the reads were aligned to the eukaryotic rRNA database available in sortmeRNA (40). The processed reads from different samples were then mapped to the BaGv2 genome, using hisat2 (41) (mapping statistics in SI Appendix, Table S6), resulting in bam files that were sorted by genomic positions, using samtools (31). They were used as inputs in StringTie (41) to create the initial transcriptome assembly with 71,042 transcripts, which was used to annotate the genome using Maker v.3 (42), resulting in 18,196 genes with 29,389 transcripts.

RNA-Seq DE Gene Analysis.

A read count matrix of the annotated genes was obtained for the samples using StringTie (41). We used the GO terms to filter out any ribosomal genes before obtaining the read counts. This approach led to the removal of 496 genes to a final set of 17,700 genes, which was used throughout the analysis. Correlations between the replicate samples were analyzed using DESeq2 (8) with a sample distance matrix. One of the antennal samples was removed due to its poor correlation with its other biological replicates. The remaining samples were used for the downstream analyses (SI Appendix, Fig. S29).

Identifying Eyespot-Specific DE Genes.

To identify eyespot-specific genes, a pairwise DE analysis was performed between eyespot and control adjacent tissues, Nes1 and Nes2, using DESeq2 (Fig. 1A and SI Appendix, Fig. S3). Common genes up-regulated and down-regulated between eyespot vs. Nes1 and eyespot vs. Nes2 with an adjusted P value (padj) of 0.05 were chosen as eyespot-specific DE genes (Dataset S1).

RNA Hierarchical Sample Clustering.

In order to identify the tissue with the closest gene expression profile to eyespots, we used all tissue samples except the eyespot control tissue samples. DE analysis between the multiple tissues was performed, using run_DE_analysis.pl script provided in Trinity tool, using DESeq2 as the method of choice for this analysis (43). Hierarchical clustering was performed for the different tissues, using genes that showed a log2fold change of |2| and padj value of 0.001, as in Fisher et al. (44). Clustering was performed using an Euclidean distance matrix derived using the DE genes for the tissues with the hclust function in R (45). The pvclust package (9) in R was used to calculate the uncertainty in the hierarchical clustering with a 1,000 bootstrap value. PCA and Spearman correlation analysis for 3-h pupal tissues were carried out using DESeq2 and corrplot (46).

ATAC-Seq Library Preparation.

We prepared ATAC libraries for the same set of tissues as we did for the RNA-seq experiment, except for the eyespot control tissues (SI Appendix, Table S7). Library preparation failed for a few groups leading to two to four biological replicates per group. Tissues were collected, flash frozen in liquid nitrogen, and stored in −80 °C, before we extracted nuclei and prepared the libraries. We used 10 to 25 individuals and ∼80,000 nuclei per replicate. Libraries were prepared as described in the Omni-ATAC protocol (47) with slight modifications. Individual tissues extracted at different time periods during the process were randomized and pooled into each replicate before extracting the nuclei. The tissues were thawed and homogenized in 2 mL of ice cold 1× homogenization buffer (HB) in a 2-mL glass douncer. Homogenization was performed by 10 strokes with pestle A, followed by 15 strokes with pestle B. The homogenized mixture was left on ice for 2 min before filtering it through a 100-µm nylon mesh filter into a DNA “low bind” 2-mL Eppendorf tube (Z666556-250EA). The filtered mixture was centrifuged at 2,500 rpm, and the pellet (the nuclei) was collected along with 50 µL of the solution at the bottom, keeping unwanted cytoplasmic RNAs in the top layers. The filtered nuclei were diluted in ATAC-resuspension buffer (RSB buffer), and 10 µL of the solution was used to count the nuclei, using a hemocytometer. Approximately 80,000 cells were used for each replicate to prepare the libraries. The tagmentation enzyme (TDE1) was obtained from Illumina (Illumina tagment dna tde1 enzyme and buffer smaller kits No. 20034197). As the concentration of the TDE1 and cell number greatly affect the identification of open-chromatin regions, we estimated the amount of enzyme needed for each reaction, using the formula: volume of enzyme = genome size of B anynana [475MB] * number of cells [80,000] *2.5/(genome size of humans [3200MB] *50,000). We used 0.65 µL (final concentration of 10.4 nM) of enzyme for each reaction. The Omni-ATAC transposition reaction was carried out as follows: 80,000 cells suspended in ATAC-RSB buffer were centrifuged at 2,500 rpm for 10 min at 4 °C. The supernatant was removed, and the nuclei-containing pellet was kept. To perform the cell lysis, 50 µL of ice-cold ATAC-RSB was added to the pellet, along with 0.1% Nonidet P-40, 0.1% Tween-20, and 0.01% digitonin. The mixture was incubated for 3 min on ice. Subsequently, 1 mL of ATAC-RSB buffer containing only 0.1% Tween and no Nonidet P-40 nor digitonin was added, and the mixture was centrifuged at 2,500 rpm. The supernatant was discarded, and the pellet was retained, to which 50 µL of transposition mixture (6.5 µL 2× TD buffer, 0.65 µL transposase (10.4 nM final concentration), 16.5 µL PBS, 0.5 µL 1% digitonin, 0.5 µL 10% Tween-20, 25.35 µL H₂O) was added. The reaction was incubated for 25 min at 37 °C at 1,000 rpm in a thermomixer. After the transpositions and tagmentation occurred, the samples were prepared for sequencing by adding Illumina/Nextera adapters with dual indexing and further PCR amplified for 10 cycles. The PCR products were purified, using a Zymo-DNA Clean & Concentrator-5 kit, and the DNA fragments were size selected between 50 and 1,500 bp, using the ProNex Size-Selective Purification System (NG2002) from Promega. The samples were sequenced, using NovosEq. 6000 with an average read depth of 30 million and 2 × 50 bp paired end reads by AIT Novogene, Singapore.

ATAC-Seq Peak Calling.

ATAC-seq analysis was perform as described in Lewis et al. (48) and Lewis and Reed (49) with modification. ATAC reads were processed, using bbduk scripts from bbmap tools to remove any adapters. The reads were mapped to the BaGv2 genome, using bowtie with the -x 1500 and -m1 parameters. Only reads with insert sizes of 1,500 bp or less and only those mapping to a unique region of the genome were mapped. Reads mapped to the mitochondrial genome were removed, using samtools idxstats, and reads marked for PCR duplicates were also removed, using GATK MarkDuplicates. We kept only paired-end mapped reads with a phred quality score of Q20 and above for downstream analysis. Because the Tn5 transposase binds to DNA as a dimer and inserts adapters of 9 bp in length at the insertion sites, the start sites of the mapped reads were adjusted to an offset of +4 bp in the forward strand and −5 bp in the reverse strand. The bam files were converted to bed files, using Bedtools (32), and we used F-SEq (50) to call peaks for each sample. Bedtools intersect was used to identify the common set of peaks for each tissue type. Peaks from all samples were merged if they were separated by 50 bp, using Bedtools merge to create 313,425 consensus peaks used for the downstream analyses. FeatureCount from the Subread package (51) was used to extract a read count matrix corresponding to the consensus peaks for all samples. The FRiP score, which is defined as the fraction of all reads that are mapped to peaks across the entire genome, was used to measure the quality of the ATAC-seq data. Our ATAC-seq data showed a median FRiP score of 0.846, which is higher than the ENCODE standard (>0.3) for the fraction of reads falling into peaks (SI Appendix, Table S8). DeepTools (52) was used to access the sample correlation between the replicates and quality of the libraries (SI Appendix, Fig. S30).

ATAC-Seq Differential Peak Analysis.

Differential peak analysis was performed using DESeq2 for 3-h pupal tissues. Peaks were considered differentially accessible with a padj value of 0.05. We also mapped the Dll319 peak identified from the FAIRE data to the BaGv2 genome, using blastn, to identify its position in the new genome assembly and test whether the ATAC-seq analysis was also able to identify it. To identify potential CREs for sal, ATAC peaks from 3-h pupal tissues were visualized using IGV and we targeted one potential candidate region (sal740) within the intronic region of sal gene loci, which is open across almost all of the tissues. Spearman correlation analysis between the 3-h pupal tissues was performed using deepTools (52).

Hi-C Analysis and Virtual 4C.

Hi-C analysis was performed as described in Lewis et al. (48).The Hi-C library used for scaffolding the B. anynana genome was reanalyzed, using the Dll319 and sal740 region as bait, to verify whether these regions interacted with the promoter of Dll and sal, respectively. Libraries were mapped to the BaGv2 assembly, using Juicer (53). We used the contact map obtained from Juicer to construct a virtual 4C plot for the window around the Dll319 and sal740 regions by placing reads in a 3-kb bin, using the script from ref. 54.

Screening and Genotyping sal740 Crispants.

Caterpillars that emerged were carefully screened under the microscope for any defects in their body, especially in the head region where sal expression is observed. Individuals showing any abnormalities were imaged and grown separately in a cup, whereas all others were grown in a separate cage. Adults were immediately frozen at −20 °C after emergence and screened later under a microscope for any defects in eyespots, wings, legs, and antennae. DNA was extracted from the crispant wing. PCR amplified the target region and amplicon sequencing was done using MisEq. 2 × 250 bp with 50,000 reads at Genewiz, China.

Hi-C Genome Assembly.

Eggs were collected from a single pair of mated B. anynana butterflies and reared. Eighteen female siblings were harvested at the wandering stage for DNA extraction. Guts were removed, and the samples were immediately flash frozen in liquid nitrogen and stored at −80 °C before the samples were sent to Dovetail Genomics to perform Chicago and Hi-C library preparation and analysis. The Chicago library preparation uses in vitro chromatin fixation, digestion, and cross-linking of regions in the genome that are close to each other in terms of three-dimensional (3D) chromatin architecture. In order to sort and scaffold the genome, 233 million reads (2 × 150 bp) were sequenced from the Chicago library and mapped to the previously published B. anynana genome (v1.2) with 10,800 scaffolds (55). The HiRise pipeline was used to identify misassemblies, to break the scaffolds, and to sort the scaffolds. Only scaffolds greater than 1 kb in length (n = 5,027) were used because scaffolds needed to be long enough for the read pairs to align and be scaffolded in accordance with the likelihood model used by HiRise. Next, 153 million reads (2*150) sequenced from the Hi-C library were mapped to the genome assembly output generated from the Chicago-HiRise pipeline to identify any misassemblies from the Chicago pipeline and correct them to produce a final genome assembly of high contiguity.

The genome assembly obtained from the HiC pipeline was ordered, using the available linkage information from Beldade et al. (56), using Chromonomer (57). A total of 289 SNP FASTA sequences were mapped to the Hi-C assembly, using blastn to identify their corresponding positions in the Hi-C genome. Using the SNP position obtained from blastn, a list describing the genetic map was manually created, which later passed through Chromonomer to sort the Hi-C assembly resulting in the final assembly (BaGv2) that was used for the current study. The BUSCO score (58) was used to check for the completeness of the gene sets in the assembly.

Genome Annotation.

The genome was repeat masked for transposable elements, small repeats, and tandem repeats before annotation as described in ref. 55. The soft repeat-masked genome was annotated, using four rounds of Maker v.3 (42). The transcriptome assembled from the RNA-seq data and gene sequences annotated from the previous version of the genome were combined and used as transcripts for the species, with transcriptome and protein sequences from Pieris rapi, J. coenia, and Bombyx mori as relative transcripts and protein homology evidences for the first round of gene predictions. Output gene predictions from each round were used as input for the next round. Snap and Augustus were used for the second round of gene predictions, followed by Genemark for the third round of gene modeling. Then we performed one final round of Snap and Augustus predictions. The minimum length of 35 amino acids was set for gene predictions. The predicted gene models were kept for genes that had an annotation edit distance (AED) score of <1 and/or had a gene ontology obtained from Interproscan (59). This resulted in 18,189 genes with 29,490 transcripts. In order to correct the annotations and produce a standardized gff3 file, the gff file obtained from Maker was run through agat_convert_sp_gxf2gxf.pl script, which is a part of AGAT tools (60). This step resulted in the removal of 82 identical isoforms and added the missing gene features, leading to a total of 18,196 genes with 29,389 transcripts. Functional annotation was performed by locally blasting the transcripts against a nonredundant (nr) protein database, using diamond blast (61), and a gene ontology analysis was performed using Interproscan in Blast2Go (62). Finally, the blast results were merged with the interproscan results in Blast2Go to produce a final functional annotation for the genome.