DOI: 10.5061/dryad.n2z34tn2x
title: Draft Papilio alphenor assembly and annotation

description: Novel phenotypes are increasingly recognized to have evolved by co-option
 of conserved genes into new developmental contexts, yet the process by
 which co-opted genes modify existing developmental programs remains
 obscure. Here we provide insight into this process by characterizing the
 role of co-opted doublesex in butterfly wing color pattern development.
 dsx is the master regulator of insect sex differentiation but was co-opted
 to control the switch between discrete non-mimetic and mimetic patterns in
 Papilio alphenor and its relatives. We found dynamic spatial and temporal
 expression pattern differences between mimetic and non-mimetic butterflies
 throughout wing development. A mimetic color pattern program is switched
 on by a pulse of dsx expression in early pupal development that
 causes acute and long-term differential gene expression, particularly in
 Wnt and Hedgehog signaling pathways. RNAi suggested opposing, novel roles
 for these pathways in mimetic pattern development. Importantly, Dsx
 co-option caused Engrailed, a key transcription factor target of Hedgehog
 signaling, to gain a novel expression domain early in pupal wing
 development that is propagated through mid-pupal development to specify
 novel mimetic patterns despite becoming decoupled from Dsx expression
 itself. Altogether, our findings provide multiple views into how co-opted
 genes can both cause and elicit changes to conserved networks and pathways
 to result in development of novel, adaptive phenotypes. This dataset
 contains the genome assembly and annotation described in the associated
 manuscript. Sequencing data, including all RNA-seq data, is available
 under the NCBI under BioProject PRJNA882073.
description: Please see our publication for further
 details. The published
 <em>Papilio polytes</em> reference genome (RefSeq
 GCF_000836215.1) was generated using <em>P. polytes
 polytes</em> from Japan, while the RNA-seq data generated in this
 study is from <em>P. alphenor</em> from the Philippines. These
 two groups diverged ~1.7 million years ago and have ~5.1% nucleotide
 divergence, resulting in low mapping rates of our
 <em>alphenor</em> RNA-seq data to the published reference. The
 <em>alphenor</em> genome assembly in this dataset is based on
 the <em>P. polytes polytes </em>assembly and a recent
 chromosome-level <em>P. bianor</em> assembly. We assembled a draft <em>alphenor
 </em>genome using PE100 sequencing from 29 individuals (BioProject
 PRJNA234541, excluding SRR1118138) and Platanus v2.0.2. Before assembly,
 we trimmed raw reads using TrimGalore!, then removed reads containing
 over-represented sequences using FastQC and an available python script. We
 then assembled all processed reads using the default platanus2 pipeline.
 Next, we assigned raw <em>alphenor</em> scaffolds to the
 RefSeq <em>polytes</em> assembly using RagTag v1.0.1, then
 assigned these <em>alphenor</em> pseudoscaffolds and all other
 <em>alphenor</em> scaffolds that hit insect sequences in the
 NCBI nr database to <em>P. bianor</em> chromosomes, again
 using RagTag v1.0.1.  Finally, we
 assembled <em>alphenor</em> <em>doublesex</em>
 alleles separately and substituted them into the final chromosome-level
 assembly. These additional steps were necessary because the individuals
 used for the full Platanus assembly were a mix of homo- and heterozygotes
 at the <em>dsx</em> locus – while the majority of the genome
 is homogeneous among these samples, the <em>dsx</em> alleles
 are very divergent except for the <em>dsx</em> coding region.
 We assembled all non-mimetic female samples BioProject PRJNA234541 using
 Platanus 2 and assigned scaffolds to the <em>polytes</em>
 assembly using RagTag. We then pulled out the <em>polytes</em>
 region corresponding to the <em>dsx</em> inversion, defined by
 Nishikawa et al. (2015) as H_locus_nonmimetic_H_scaffold:1931762-2054949,
 and used it to replace the corresponding region in the chromosome-level
 <em>alphenor</em> assembly. Similarly, we assembled the
 mimetic <em>dsx</em> allele by assembling all mimetic female
 samples from BioProject PRJNA234541 using Platanus 2, assigning those
 scaffolds to the RefSeq assembly, and extracting the
 H_locus_mimetic_hetero scaffold. We added this scaffold as chrH in the
 final <em>alphenor</em> assembly.  We assembled the <em>alphenor</em>
 mitochondrial genome using NOVOplasty v4.2 using sequencing data from
 SRR1108726 and the RefSeq mtDNA assembly for polytes (NC_024742.1) as the
 seed sequence. This resulted in a single circularized sequence of 15,247
 bp. We annotated the
 <em>alphenor</em> genome using EvidenceModeler 1.1.1. We first
 assembled a high-quality transcript database using PASA, our SE50 data,
 and PE100 and SE50 data from Nallu et al. (2018). After adapter trimming,
 we performed de novo and genome-guided assembly using Trinity v2.10.0 and
 genome-guided assembly using StringTie v1.3.3. RNA-seq data was also
 mapped to <em>alphenor</em> chromosomes using STAR 2.6.1d, and
 the resulting alignments were used to generate genome-guided assemblies
 with Trinity and StringTie 1.3.1. We combined de novo and genome-guided
 assemblies using PASA 2.4.1. Evidence for protein-coding regions came from
 mapping the UniProt/Swiss-Prot (2020_06) database and all Papilionoidea
 proteins available in NCBI’s GenBank nr protein database (downloaded
 6/2020) using exonerate. We identified high-quality multi-exon
 protein-coding PASA transcripts using TransDecoder
 (transdecoder.github.io), then used these models to train and run
 Genemark-ET 4 and GlimmerHMM 3.0.4. We also predicted gene models using
 Augustus 3.3.2, the supplied heliconius_melpomene1 parameter set, and
 hints derived from RNA-seq and protein mapping above. Augustus predictions
 with &gt;90% of their length covered by hints were considered
 high-quality models. Transcript, protein, and ab initio data were
 integrated using EVM.  Raw EVM
 models were then updated twice using PASA to add UTRs and identify
 alternative transcripts. Gene models derived from transposable element
 proteins were identified using BLASTp and removed from the annotation set.
 The final annotation comprises 17,342 genes encoding 26,991 protein-coding
 transcripts, containing 95% and missing 3% of endopterygota single-copy
 orthologs according to BUSCO v3 and OrthoDB v9. We functionally annotated
 protein models using eggNOG’s emapper-2.0.1b utility and the v2.0 eggNOG
 database. Transcript and protein sequences were
 extracted from this final annotation file using Cufflinks' gffread
 utility.
description: Standard Unix utilities such as less, grep, and emacs can be used
 to visualize these text files. Fasta sequences can be retrieved using
 samtools' faidx utility. Genome fasta and GFF3 annotation files can
 be directly loaded into a genome viewer such as IGV.