Positional cloning

All Limenitis specimens utilized for genetic linkage mapping were collected from a single locality in Pennsylvania (Supplementary Table 2). Wild-captured, mated female specimens were captured and fed a mixture of honey and water twice daily, and were secured on Prunus serotina and/or Salyx babylonica to encourage oviposition. Larvae were raised directly on host plants with one brood per enclosure. Pupae were collected and placed in labelled containers to prevent adults from mating upon emergence. Adult butterflies were transferred to envelopes numbered with their sibling group and according to their order of emergence. Adults were then photographed and crossed via hand-pairing41. Mapping families were generated via backcrossing heterozygous mimetic males to fully banded, non-mimetic females (homozygous recessive at the major gene controlling mimicry), as linkage maps constructed with heterozygous females are uninformative because there is no recombination during oogenesis42. Following mating, the wings and tissues of male butterflies were immediately archived, and the wings and tissues of female butterflies were archived when oviposition was ceased. Progeny from mapping crosses were photographed upon emergence, and wings and tissues were archived.

Medial banding in Limenitis is controlled by two, incompletely dominant, alleles at a single locus43, and at least one dominant modifier that influences the penetrance of the white-banded allele in heterozygous individuals. All progeny from mapping crosses displayed either the mimetic phenotype (heterozygous with dominant modifier) or the white-banded (homozygous) phenotype. We scored wing patterns of the resulting 111 progeny of the mapping family based on the presence or absence of hindwing and forewing medial white banding. Forty-five were mimetic (23 males and 22 females) and 66 were fully banded (30 males and 36 females; χ2=3.973, two-tailed P value=0.046). Although these values differ weakly from our statistical expectations, the results of numerous other crosses carried out by SPM over the last 10 years support a model of Mendelian inheritance.

Wings were removed from each butterfly, photographed and archived in glassine envelopes. Remaining whole bodies were placed inside 1.5 ml microcentrifuge tubes with 100% ethanol to preserve the genomic DNA. Wing muscle tissue was dissected from each archived butterflies, and DNA extractions were performed using the DNeasy Blood and Tissue Kit (Qiagen, Inc.). Amplified fragment length polymorphism (AFLP) genotypes were generated using the AFLP Plant Mapping kit for small plant genomes (Applied Biosystems, Inc.) following the manufacturer’s instructions. In brief, fragments were generated using restriction enzymes EcoRI and MseI. Next, fragments were ligated to adaptors and selectively amplified during two separate rounds of PCR. Once the PCRs was completed, reactions for each fragment were spiked with an internal lane size standard to ensure proper size matching of AFLP across samples. Fragment analysis was performed on an 3730xl DNA Analyzer (Applied BioSystems, Inc.). AFLP electropherograms were analysed in Applied Biosystem’s GeneMapper software, version 3.7 (Applied Biosystems, Inc.), with AFLP alleles scored as either present or absent and according to fragment size. Peak height threshold was set to 80 RFU. Fragments smaller than 50 base pairs and larger than 550 base pairs were excluded from analysis. The parents and progeny (n=111) were then genotyped for 60 unique AFLP primer pair combinations44 analysed separately in GeneMapper, version 3.7 (Applied Biosystems, Inc.), and any marker that had a large number of missing genotypes was removed from the data set. Of the resulting 2,571 AFLP, 506 (19.7%) were diagnostic, mappable markers (present in the male, absent in the female and segregating in the progeny in a near 1:1 pattern).

AFLP marker order and position were estimated using JoinMap 3.0 (ref. 45), which used χ2-tests for segregation distortion and tests of independence to filter spurious genotypes. Log likelihood ratio scores above 3.0 for markers are considered significant evidence for linkage (95% probability). Markers were separated into linkage groups and are selected for mapping at a minimum log likelihood ratio score of 4.0. Final map distances were corrected using Kosambi’s mapping function; this function is designed to account for the observation that larger chromosomes are more likely to have double crossovers than small chromosomes and results in shorter, more accurate linkage maps compared with those maps using Haldane’s mapping function that assumes complete interference45. Mapping of the diagnostic markers resulted in 30 linkage groups, varying in length from 46 to 107 cM, and in agreement with chromosomal number estimates for this species46. The total map length was 2,250 cM with an average distance of 4.5 cM between AFLP markers. Given the genome size estimate for Limenitis47 of 388 Mb (+/− 7 Mb for females and 4 Mb for males), each cM in our map represents ~160.4 kb.

Direct comparisons among linkage maps based on dominant markers, such as AFLPs, are not possible. Therefore, we generated diagnostic SNP markers from highly conserved nuclear genes by designing Limenitis-specific primers from annotated expressed sequence tags (ESTs) developed via 454 pyrosequencing of 48 h pupal wing disc cDNA, and genotyped the entire backcross mapping brood for each SNP by Sanger sequencing to facilitate comparisons with other published Lepidopteran genetic maps, with an emphasis on Heliconius. Such comparisons are greatly facilitated in Lepidoptera because macrosynteny is highly conserved33,48,49, thereby allowing direct identification of homologous chromosomal linkage groups across systems.

The final linkage map contains 57 SNPs from 54 conserved nuclear genes, such as 506 AFLPs, 1 BAC end sequence and 2 colour pattern loci (wbd and ird), providing 565 distinct markers with a total map distance of 2,248 cM across 30 linkage groups (only the W chromosome was not mapped). The average linkage group was 74.9 cM. Both the total map length and the average linkage group length remained nearly the identical to the dominant marker map. However, the addition of 57 new SNPs reduced the average distance between markers from 4.44 to 3.99 cM, an estimated 6,000 bp difference per cM. Importantly, mapping of the nuclear genes demonstrated that mimetic variation in Limenitis segregates with several genes known to be linked to the chromosome in Heliconius, which is known to house the colour-patterning locus, Ac50. The addition of 19 syntenic nuclear markers49 known to be linked to Ac in Heliconius (Supplementary Fig. 1 and Supplementary Table 1) resulted in a more accurate positioning of markers on the Limenitis wing patterning chromosome, and reduced the zero-recombinant mapping interval to a 2-cM region containing two adjacent chitin synthase genes and the candidate gene, WntA20.

Developmental patterns of gene expression and functional tests

Mimetic and non-mimetic female Limenitis were wild captured from Baltimore county, Maryland, and White Mountain National Forest, New Hampshire, respectively. Lab-reared offspring from true-breeding individuals were mated to full-sibs, and the resulting progeny were utilized as reported for all following developmental and functional tests.

Pupae aged 8–16 h after pupation were injected with 10 or 20 μg μl−1 heparin or sterile H 2 O using a pulled glass micropipette mounted on a 10-μl cut pipette tip and a 2–20- μl pipette. Pupae were surface sterilized with ethanol and injected on the left side in an interstice that separates the baso-posterior parts of the developing forewing and hindwing. As in previous reports20, H 2 O controls showed no effects, heparin had systemic effects on both left and right wing patterns, and control and heparin injections did not produce local damage artefacts. All of the results presented in Fig. 2 were replicated at least three times per morph and dosage (including H 2 O control).

Individual mimetic and non-mimetic Limenitis were sampled at the 5th instar stage of development (Supplementary Fig. 2A). Each individual was anaesthetized in ice-cold water and their wing discs dissected in PBS, incubated in cold fixative (formaldehyde 9% in PBS containing 50 mM EGTA) for 30–35 min, rinsed in PBS with 0.01% Tween20 (PBST), and then dehydrated with increasing concentrations of methanol and kept for long-term storage in methanol at −20 °C. For in situ hybridization, we followed the method described in ref. 20. In brief, wing discs preserved in methanol were rehydrated in increasing concentrations of cold PBST, washed in cold PBST and freed from their peripodial membrane using fine forceps. The tissues were then post fixed 20 min on ice in PBS containing 5.5% formaldehyde, transferred to a standard hybridization buffer, incubated in supplemented hybridization buffer for 16–40 h at 62 °C. Tissues were washed eight times, and, then for secondary detection of the riboprobe, the tissues were blocked and then incubated with a 1:4,000 dilution of anti-digoxigenin alkaline phosphatase Fab fragments (Roche Applied Science, Indianapolis, Indiana, USA). Tissues were again washed 10 times for 10–120 min, and finally stained with BM Purple (Roche Applied Science) for 4–8 h at room temperature. Stained tissues were then washed in PBST 2 mM EDTA and slide mounted in PBS containing 60% glycerol. mRNA in situ hybridizations were photographed with a Nikon Coolpix P5100 digital camera (Nikon Inc., Melville, New York, USA) mounted with a LNS- 30D/P51 adapter (Zarf Enterprises, Spokane, Washington, USA) on a Leica S4E microscope (Leica Microsystems, Buffalo Grove, Illinois, USA).

Mimetic and non-mimetic individuals were captured in Massachusetts and New Hampshire, respectively, and allowed to oviposit in the laboratory. Mimetic and non-mimetic individuals were then crossed to siblings to ensure that they were true-breeding for each respective phenotype. Individuals from each cross were sampled at the 5th instar (n=3 mimetic and n=3 non-mimetic), prepupa (identified by stereotypic ‘j-curve’ hanging from leaves; n=3 mimetic and n=3 non-mimetic), early pupation (<48 h post pupation; n=3 mimetic and n=3 non-mimetic), and late pupation (>48 h post pupation; n=3 mimetic and n=3 non-mimetic). Wing discs were dissected from each individual in cold PBS, wing discs were then stored in RNAlater (Ambion, Inc.) following manufacturer’s instructions. RNA was extracted using an RNAeasy-kit (Qiagen, Inc.), and individual RNAseq libraries were prepared using the TrueSeq RNA sample prep kit (Illumina, Inc.) at the Michigan State University Genomics core facility. Each library pooled, and pools were sequenced across two lanes of an Illumina HiSeq 2500, using 2 × 150 bp reads.

Individual reads were quality filtered and trimmed using custom python scripts, then aligned to BAC sequences representing mimetic and non-mimetic haplotypes using the TopHat pipeline51. TopHat identified genes and exons automatically based on read alignment. Differential expression of whole genes was determined using Cufflinks51, and differential expression of exons was determined using DEXSeq30. Differential expression of the WntA was not detected at any developmental stage. Differential expression of WntA exon 1 was detected between mimetic and non-mimetic individuals only in the 5th instar stage, where WntA expression is highest (P<0.0001 after Benjamini–Hochberg correction for false positives).

Generation of BAC tile path

A custom Limenitis BAC library was constructed from eight individuals captured from a single hybrid population from Pennsylvania (Supplementary Table 1). Nylon arrays of BAC clones were prepared by Amplicon Express (Pullman, WA), and were screened using Mat1 and WntA probes. Of the entire library of ~18,000 clones, 10 BACs screened positive for Mat1, 9 BACs screened positive for WntA and 1 clone, 21G6, screened positive for both markers. These 19 clones were selected for end sequencing and high-depth next-generation Illumina sequencing and assembly, performed by Amplicon Express. Because of considerable structural variation between these BAC contigs, sequences were aligned in multiple steps. First, BAC end sequences of all 19 clones were mapped using BLAST to each individual contig, and clones were aligned manually. Next, when a tile of multiple contigs spanning the zero-recombinant interval was determined, we performed a final alignment of the constituent BAC contigs using the MAFFT52 algorithm as implemented by Genious Pro (Biomatters, Inc.). The final reference was comprised of three major contigs representing clones from 21G6, 60N10 and 70O17. A fourth BAC contig, 43D10, contained the alternative Limenitis haplotype reported in the text. The assembled reference sequence for the zero-recombinant interval (‘Limenitis AC scaffold’) used for our comparisons contained the haplotype typical of the mimetic, unbanded phenotype. A FASTA file containing all BAC sequences is available as a supplementary file (Supplementary Data 4).

Population genomic resequencing

We collected adult butterflies for DNA sampling from a variety of locations (Supplementary Table 1), representing 120 individual Limenitis spp from sites in Georgia, Virginia, Vermont, Pennsylvania and Maine, as well as 25 H. c. alithea from sites in Ecuador, and finally 10 additional localities in Costa Rica, representing H. c. galanthus (n=10) and H. pachinus (n=10). All samples were euthanized, and then wings were separated and placed in glassine envelopes while the bodies were stored in 95% ethanol. For each sample, we extracted DNA from butterfly wing muscles using DNAeasy kits (Qiagen, Inc.) following manufacturer’s instructions.

For Limenitis and H. c. alithea, we performed short-read genomic library construction using the Nextera Sample Preparation kit, following the manufacturer’s instructions. These libraries were sequenced using an Illumina HiSeq 2000, with 2 × 150 bp paired-end reads, at the Bauer Sequencing Core Facility (Harvard University) to ~10 × coverage per individual. We subsequently generated ~8 × coverage for both parents of our backcross mapping brood using 2 × 250 bp MiSeq runs. For the remaining Heliconius spp., a custom Illumina sequencing library with a 500-bp insert was prepared for each sample and sequenced to an average depth of 16 × coverage using an Illumina HiSeq 2000 (2 × 100 paired-end sequencing). Library preparation and sequencing were performed at the Beijing Genome Institute. This sequencing data is archived in the NCBI Short Read Archive with the following BioProject accession number: PRJNA226620. All remaining Limenitis and H. c. sequencing data is archived in the NCBI Short Read Archive with the following accession BioProject accession number: PRJNA252628.

BAC sequence alignment

Limenitis spp: Following sequencing, sequence libraries representing each sample were trimmed and filtered for quality using the fastx toolkit (CSHL), scythe and sickle (UC-Davis). Next, filtered reads were aligned to the Limenitis reference BAC tile path (see the section ‘Generation of BAC tile path’) using the ‘very fast’ parameter set with local alignment. SNPs were called simultaneously for all samples using the multi-allelic calling function in GATK version 1.5 (refs 53, 54). Positions with a total SNP quality <40 were filtered from subsequent analyses. This resulted in a data set containing 9896 SNPs. A custom perl pipeline was used to identify biallelic SNPs within this data set that were fixed between phenotypes.

Average pairwise LD, for SNPs was calculated as the squared correlation coefficient (r2) between allele counts observed between each SNP and its 200 nearest neighbours (5,000 bp, maximum distance) using the VCFtools software package55. Average r2 at this location was then was then calculated for all r2 values, and the process was repeated. This approach is computationally feasible for large data sets since it does not require haplotype reconstruction, but it provides only an approximation of the true LD56. Based on these analyses, we identified a large (30 kb) segregating haplotype that was perfectly association with wing pattern phenotypes in Limenitis (Supplementary Fig. 3).

Hybrid zone transect genotyping

To further investigate the phenotypic association between Limenitis wing pattern morphs and the segregating haplotype, we designed a TaqMan probe-based gene expression assay using custom primers. The two probes, which were haplotype specific, were labelled with a FAM or VIC dye label on the 5′ end, and a minor grove binder non-fluorescent quencher on the 3′ end. We PCR amplified all individuals using a 20-μl volume, following the manufacturer’s suggested conditions, and then took end point fluorescence measurements, using an Eppendorf RealPlex2 mastercycler, to call heterozygotes and both homozygote genotypes. The frequency of each haplotype relative to geographic sampling locality is shown in Fig. 4. The presence of heterozygous individuals in populations outside of the hybrid zone reflects geographic differences in the frequency of alleles at the modifying locus. Previous crossing experiments indicate that the dominant modifier is at high frequency near the hybrid zone but declines in frequency with increasing geographic distance from the hybrid zone.

Population structure tests

To test the null hypothesis of geographic population structure, non-mimetic (L. a. arthemis) and mimetic (L. a. astyanax) butterflies (n=417 individuals) were sampled across two independent transects of the hybrid zone (n=10 populations/transect; Supplementary Table 3), and genotyped with 12 selective AFLP primer combinations following the same protocol described above (see Genotyping). The 12 unique AFLP primer pair combinations resulted in the presence or absence of 2,723 AFLP loci, of which 490 AFLPs had a minor allele frequency ≥0.10 and were retained for subsequent analyses (Supplementary Table 4). Outlier analysis was then performed to identify loci experiencing selection using the programme BAYESCAN v. 2.0 (ref. 57). Within BAYESCAN, estimation of model parameters was tuned automatically on the basis of short pilot runs (10 pilot runs, length 5,000), using default chain parameters. Loci were then ranked according to their estimated posterior probability, and all loci showing log(Bayes Factor)>2 ((P[ai] 6¼ 0)>0.99) were treated as outliers (Supplementary Fig. 5). Partitioning the AFLP data into outlier (n=62; FDR=0.001) and non-outlier (n=428) loci revealed low overall population differentiation (F ST neutral=0.09 versus 0.51 for outliers) and no evidence for population structure based on wing pattern. We also used the programme STRUCTURE v. 2.2 (ref. 58) to cluster individuals from each transect on the basis of their multilocus AFLP genotypes. For both transects, clustering runs (burn-in=5,000 repetitions and main run=5,000,000 repetitions) using the admixture model, and allowing the number of clusters to vary from K=1 to 10, returned the maximum log likelihood value for K=4. However, inspection of the clustering results indicates no correspondence between individual genotypes and geography or wing pattern (Supplementary Fig. 6).

Ac reference sequence alignment

H. c. alithea: The described procedure for Limenitis was followed, except sequences were aligned to the H. melpomene scaffold known to contain the WntA gene20. The NCBI GenBank accession number for this scaffold is HE668478. This resulted in a data set containing 37,347 SNPs.

Additional Heliconius species: Reads were trimmed and filtered for quality as described above. These additional Heliconius species were subsequently aligned to the Hmel 1.1 reference genome33 using Stampy55. SNPs were called simultaneously for all Heliconius spp samples using the multi-allelic calling function in GATK version 1.5 (refs 53, 54). Positions with a total SNP quality <40 were filtered from subsequent analyses. This resulted in a data set containing 42,522 SNPs. We found a single 1.8-kb indel that was perfectly associated with allelic variation in Ac phenotypes among H. c. galanthus, H. pachinus and the two phenotypes of H. c. alithea (Supplementary Fig. 8).