dataesr/scientific-paragraphs-categorization huggingface dataset:
Combined with already available resources, which we also improve both in terms of assembly contiguity and gene annotation, we generate a genomic dataset that comprises ~75% of all the species in the Heliconiini tribe, to our knowledge one of the most comprehensive efforts to sample an insect tribe at high taxonomic density. With this, we produce a comprehensive dated phylogeny for Heliconiini and explore patterns of gene flow across the tribe. We test if a significant and substantial amount of genomic change occurred not only at the stem of Heliconius, but also at more basal branches within the Heliconiini tribe, pre-dating the range of innovations seen in Heliconius. Finally, we investigate structural and adaptive aspects of genome evolution across the radiation and during key ecological transitions, and explore evidence of accelerated evolution in putative regulatory elements.
Improved resolution of phylogenetic relationships and signatures of introgression across the genome To generate the species tree, we first compiled a total data set of 4,011,390 base pairs of aligned protein-coding DNA obtained from the single-copy orthologous groups (scOGs) obtained from a high quality and revised genomic resources (Supplementary Note 1, Supplementary Figs. 1-19). The alignment has over 1. 5 M parsimony-informative, ~500k singleton sites, and 1. 9k constant sites. A species-level phylogeny was determined with a maximum-likelihood (ML) analysis (Supplementary Figs. 19 and 20), and used to estimate divergence dates (Fig. 1c. Supplementary Fig. 21 , and Supplementary Data 3). Although the topology is widely consistent with previously inferred phylogenetic relationships  8  , we identify differences within some Heliconius clades, where the Silvaniform and Melpomene clades are now paraphyletic, and among other genera of Heliconiini, with Podothricha telesiphe and Dryas iulia now sister lineages, outgrouped by Dryadula phaetusa and Philaethria dido.
The fraction of the genome that introgressed (average f-branch statistic) across all triplet comparisons and coalescent units (CUs), as a proxy of ILS, from each chromosome and the Z chromosome versus all autosomes (see Supplementary Note 3 and Supplementary Figs. 25, 26) show that, indeed, the Z chromosome has a lesser degree of introgression and ILS overall, but this effect does not change the topologies' frequencies in favor of Topology 3, which stays ~5% of the entire chromosome, whereas Topology 1 increase to 56%. Overall, given the methods currently available for large phylogenomic datasets such as ours (but see Thawornwattana et al. 17  ), the landscape of local history seems to confirm the species tree as the most consistent topology, with H. aoede clustering within Heliconius clades.
Evolution of genome size and content Variation in genome size can be formalized in an "accordion" model  18  where genomes gain, lose, or maintain its size in equilibrium in each species, due to a balance between expansions in transposable elements (TEs) and large segmental deletions. By reconstructing ancestral genome sizes at each node in the Heliconiini phylogeny, we found that the MRCA of Heliconiinae experienced a 30% contraction from ~406 Mb at stem of Heliconiinae to ~282 Mb for Melpomene/Silvaniform clade; while, at the same time, other branches leading to Philaethria, Dryadula, Dryas and Podothricha, and the Erato, Doris and Wallacei clades within the genus Heliconius, had independent expansions. Strikingly, H. aoede shows a loss of about 68 Mb, a fifth of its genome size (~22%) from its ancestral node (Fig.
Candidate genes for derived traits of Heliconius Within the Heliconiinae, Heliconius display a number of divergent traits and innovations  10. Here, we highlight how our results reveal new biological insights into these traits, focusing on two case studies; changes in neural composition in Heliconiini, and the enzymatic processes associated with breaking down pollen walls to aid their digestion during pollen feeding. These two examples, illustrate the potential of large, densely sampled genomic datasets to both generate and test adaptive gene-phenotype hypotheses, using both unguided and more targeted analyses. Within in central brain, mushroom bodies are paired organs that receive visual and/or olfactory information, and play a pivotal role in learning and memory  54. These structures show huge variation across Heliconiini, but a particularly large expansion occurred at the Heliconius stem, where mushroom body volume and neuron number more increased by several-fold, accompanied by a major shift towards increased dedication to processing visual information  11, 26.
Moe drives cortical remodeling of dividing neuroblasts  57  , while Fbxl7 affects Hippo signaling pathway activity  58. Finally, Ctr9, dachsous (ds), falafel (flfl), and locomotion defects (loco) were identified at the Heliconius stem (Supplementary Data 9). Ctr9 is involved in the proliferation and differentiation of the central nervous system  59  , Flfl is required for asymmetric division of neuroblasts, cell polarity and neurogenesis in mushroom bodies  60  , Ds is a cadherin that interacts with the Hippo signaling pathway  61  , and loco is an activator of glial cell fate, essential cells in efficiently operating nervous systems  62. Similarly, our analysis of conserved non-coding elements reveals multiple loci nearby genes with known roles in neural development, synaptic pruning, and longterm memory.
We identified 233 cocoonase loci (Supplementary Data 15) across all Heliconiinae and found that the duplications not only predate the split between Heliconius and Eueides, but affect the whole Heliconiini tribe and its outgroup S. mormonia (Fig. 7a ). All species have at least four copies, located at the minus strand of chromosome 15 with remarkable conserved synteny (Fig. 7b  and Supplementary Data 16). Substrate, cleavage, and active sites of the functional domain show very high conservation throughout the dataset. Three independent tandem duplications from the same original copy are very likely responsible of the emergence of Coc1A, Coc1B, Coc2 and Coc3 (Fig. 7a ). High level of purifying selection is detected across the four OGs (Fig. 7a ). A scan of all internal branches for signs of diversifying positive selection shows that the branches of D.
In conclusion, we have curated available genomic data and reference genomes to build a tribe-wide dataset for Heliconiini butterflies. Using the resulting phylogenetic framework, we examined patterns of genomic change at points in the species tree around which key phenotypic innovations are expected. We investigated the evolution of genome size, its effect on protein-coding gene expansions and contraction, and selective forces such as diversifying positive selection on protein coding genes and the acceleration of conserved non-coding genes. Supported by the characterization of all these genomic features, our analyses ultimately allowed us to narrow down candidate genes that could be further tested to explore more in depth the molecular architecture of key innovations in this enigmatic group of butterflies.
Finally, our test for acceleration of putative cis-regulatory elements (CNEEs) at the stem of Heliconius identified more prevalent positive selection on non-coding elements compared to protein coding genes at the origin of Heliconius. This suggests the suite of derived phenotypes in this genus might have largely evolved through changes in gene expression via modification of regulatory elements (e. g. : promoters, enhancers, and silencers)  41, 66. In conclusion, our work offers a comprehensive view to the evolutionary history of an enigmatic tribe of butterflies, the evolution of their genomic architectures, and provides the most thorough analysis of potential molecular changes linked to the physiological and behavioral innovations of a diverse group of butterflies. These gene-phenotype hypothesis, alongside our comprehensive dataset, provide new opportunities to test and derive causative links between molecular and trait innovations.
All the extractions were finally pooled keeping the same final RNA concentration from all samples. Due to different DNA quality, different 3 rd generation sequencing technologies were applied. Pacific Bioscience Sequel II / HiFi was used to sequence the S. mormonia washingtonia genome at the Genomics Resource Center (GRC), University of Maryland Baltimore (UMD-IGS). Other Pacific Bioscience data were generated from Dr. phaetusa, Di. juno, and A. vanillae (Florida) at the Center for Genomic Research, University of Liverpool using PacBio sequel SMRT cell (2. 0 chemistry), adding a library of Illumina DNAseq (HiSeq2500 150 × 2) data for D. phaetusa. For Ph. dido, Po. telesiphe, A. vanilla (Costa Rica), H. aoede, and H. doris, the 10x Chromium Library Prep was adopted and Illumina sequencing using 150 bp paired-end reads with NovaSeq FC S2, generating ~40 Gbp per species, and adding Nanopore 1D (up to 5.
PacBio genome assembly PacBio Hifi CCS reads were directly used to assemble the genome of Speyeria mormonia washingtonia (Washington, USA) using HIFIASM V0. 12-r304 [settings: -f0] 67 ; while for other species PacBio reads were corrected, trimmed and assembled using CANU V1. 8 + 356 changes  68  [settings: genomeSize = 400 m; corMhapSensitivity = normal; corOut Coverage = 100; correctedErrorRate = 0. 105; ovlMerThreshold = 500; batOptions = -dg 3 -db 3 -dr 1 -ca 500 -cp 50], as in Cicconardi et al. (2021)  20. The resulting raw assemblies were subsequently corrected by remapping all uncorrected raw PacBio reads with PBMM2 V1. 0. 0 and ARROW V2. 3. 3 (https://github. com/PacificBiosciences/GenomicConsen sus) was used to correct the assembly with three iterations. To further error correct contigs, short Illumina reads were mapped with STAR V2.
7. 2c 69 , and PILON V1. 23  70  was run for five iterations. Assemblies were then processed with PURGE HAPLOTIGS V20191008 71 using ad hoc -a parameter to optimize the removal of haplocontigs, artificially duplicated genomic regions due to heterozygosity. To correct for mis-assemblies POLAR STAR (https://github. com/phasegenomics/polar_star) was employed. This tool calculates read depth of aligned PacBio reads to the assembly at each base, smoothed in a 100 bp sliding window, merging regions of high, low, and normal depth. Low read depth outliers are identified, and contigs are broken at each such location. Contigs were then rescaffolded using P_RNA_SCAFFOLDER  72  , which uses information from RNAseq mapping (if RNAseq were available), and LRSCAF V1. 1. 5 73 , which uses information from long-reads.
Curation of available illumina assemblies We took advantage of the available Heliconiini "Discovar" assemblies  7  , including them in our dataset, but before proceeding with downstream analyses (e. g. , gene annotation), a small, but effective curation was applied. At first, we checked for contaminants (see below), and then checked for completeness (using BUSCO, see above), the presence of haplocontigs. Raw Illumina reads were remapped onto their own assembly and PURGE HAPLOTIGS, with ad hoc -a parameter, was adopted to remove them, followed by a scaffolding procedure with SOAPDENOVO2 (127mer), and further synteny mapping implemented with ALLMAPS, using their closest available relative with a reference assembly. This procedure was adopted to maximize the contiguity. For future studies it should be noted that this may inhibit further insight into different genomic rearrangements from the reference species by assuming conserved synteny, a procedure similar to this was recently adopted to scaffolded draft genomes  95.
Bacterial contamination & assembly completeness assessment After the genome assembly stage all datasets were analyzed to remove contaminants. BLOBTOOLS V1. 1. 1 96 was implemented using BLASTN [-evalue 1e-25 -max_target_seqs 1] and the NCBI nucleotide collection (#seqs: 49,266,009, retrieved September 2018). Taxon identities for each hit were retrieved and we filtered out any scaffold and contig that were assigned to fungal and bacterial contaminants. Furthermore, mitochondrial sequences were identified by blasting (BLASTN) contigs and scaffolds against the mitochondrial genome of available Heliconius ssp. (HE579083. 1, NC_026564. 1, NC_026463. 1, KP281778. 1, KP294327. 1, KP100653. 1, NC_024744. 1, NC_024741. 1, KM068091. 1, KM014809. 1, NC_027516. 1, KP784455. 1, NC_024864. 1, KM208636. 1).
We calculated rates for transposition (d P ), insertion (d I ), deletion (d D ), inversion (d V ), and duplication (d U ) events per million years (Ma) of evolution, based on the inferred divergence estimates from the dated phylogeny (see below). We subsequently calculated the ratio between d D /d U , in order to better understand genome changes over time. The evolution of genome size was also assessed using an ancestral state reconstruction method implemented in the R package PHYTOOLS 102 , using the dated phylogeny (see below). The maximum likelihood ancestral character state at each node was inferred across the phylogeny and the function CONTMAP was used to plot these continuous character traits onto the phylogeny in PHYTOOLS. To explore evolutionary conservation at individual alignment sites, PHYLOP scores were computed.
A neutral model was first generated as implemented in HALPHYLOPTRAIN. PY script (CACTUS package) using protein-coding exon coordinates of H. melpomene. Using that neutral model, a PhyloP score measure (CONACC) was computed for each species in the phylogeny. A non-overlapped sliding window of 10 bp was adopted and data partitioned according to coding, intronic, 5'UTR Gene prediction and transcriptome annotation NCBI SRA archive was explored to look for available RNA-seq data. The best SRA archives were selected, downloaded (Supplementary Data 2) and added to the data generated in this study to re-annotate genes for available annotations (H. erato v. 1 and H. melpomene v. 2. 5), and annotate genes for also for available assemblies. To do so the following pipelines was implemented (Supplementary Fig.
30 105 and AUGUSTUS V3. 4. 0 106 , along with the masked genomes, generated with REPEATMASKER V4. 1. 1 107 , using the Lepidoptera database, and the alignment of protein sequences from model species (Drosophila melanogaster, Bombyx mori, Bicyclus anynana, Danaus plexippus, H. erato, and H. melpomene) using EXONERATE V2. 46. 2  97. BRAKER is an automated pipeline that aims to predict genes using iterative training of AUGUS-TUS. GENEMARK-predicted genes are filtered and provided for AUGUS-TUS training, followed by AUGUSTUS prediction, integrating the RNA-Seq and protein alignment information, to generate the final gene predictions. To the regular BRAKER run we added a further UTR annotation step [braker. pl --addUTR=on], which uses the RNA-Seq coverage data via the tool GUSHR V1. 1. 0 to add UTRs to the augustus.
hints. gtf file; it does not perform training of AUGUSTUS or gene prediction with AUGUSTUS and UTR parameters. Unfortunately, because the generated output (augustus. hints. utr. gff3) presents many syntax errors, we tried to mitigate them with a custom python script (BRAKER-GUSHRGFF3TOBED12V1. 0. PY) available at https://github. com/ francicco/-ComparativeGenomicsOfHeliconiini. To generate the de novo transcriptome assemblies, quality filtered reads were assembled using TRINITY V2. 10. 0 108,109 separately for each tissue. The generated contigs were subsequently aligned to the genome using MINIMAP2. Coordinates for the aligned contigs were used to extract nt sequences, and TRANSDECODER V5. 5. 0 (http://transdecoder. github. io/) (minimum amino acid length >50), using homologs from UNIPROt database 110 and Lepidoptera proteome (see below) found with DELTABLAST V.
PY) available at https:// github. com/francicco/-ComparativeGenomicsOfHeliconiini. The whole procedure generated four annotations per species, a predicted one, a de novo, and two ab initio annotations. These were combined and used with STAR to map pooled reads again to their own assembly. The BAM file that was subsequently obtained was used as input for PORTCULLIS V1. 1. 2 118 [--threshold 0. 5] to generate a splice-site DB. Finally, all these elements (transcript and spice-site annotations) were combined together using MIKADO V2. 3. 3 119 [--scoring insects. yaml -bt UNIPROTDB+Lepidoptera -mode split]. For species without their own transcriptome data (RNA-seq) the procedure to annotate their genomes was as follows. For each one, all mRNA transcripts from all of the annotated species so far were mapped to its genome using MINIMAP2 and the resulting BAM file was given to BRAKER [--addUTR=on] as input.
Intron evolution analyses Intron length evolution was studied looking at both the interaction with transposable elements (TEs), and looking at the rate of gains and losses across the phylogeny. Intronic regions were extracted from the longest transcript of each gene model using the annotations, as in Cicconardi et al. 2021  20  , and their sequences were used to annotate TEs with  REPEATMASKER. The correlation between genome size and median intron length was also computed with the PGLS function in R from the CAPER package 121. The function implements GLS models accounting for phylogeny. We included a variance covariance array representing the phylogeny within the comparative dataset of dimension three (vcv = TRUE, vcv. dim = 3). The branch lengths were transformed optimizing between bounds using maximum likelihood (lambda = 'ML').
It also has the ability to detect chimeric proteins resulting from gene-fusion events and assigns these proteins to the corresponding OGs. Once parameters were optimized [-kmer_size 300 -e_value 1e-40 -sp_overlap 0. 1 -min_nb_hits 20 -min_weight 0. 3] (Romain Derelle, personal communication) a first iteration of the analysis was run using the whole dataset as input. Chimeric transcripts where then identified. For the E. isabella transcriptome the chimeras where manually curated by blasting single proteins to the BLAST server and then used to set up a custom automatic parser in python (BROCCOLICHIMERASPLITDATAGATHER. PY; available at https://github. com/francicco/-ComparativeGenomicsOfHeliconiini), before a second iteration was executed with the generated data. From this orthology table, pseudo-paralog loci from the highly fragmented genomes were identified and excluded if the sequence was shorter than the average value, with one and three standard deviations used as thresholds for scaffold and sub-chromosome level assembly, respectively.
From each OG a putative functional annotation was performed by identifying both the protein domain architecture and GO terms. For each sequence, HMMER V3. 3. 2 (HMMSCAN)  125  was first used to predict PFAM domains, then DAMA V2 (Domain Annotation by a Multi-objective Approach)  126  was applied to identify architectures, combining scores of domain matches, previously observed multidomain co-occurrence, and domain overlapping. Annotation of GO terms was assigned with two strategies: homology-and predictive-based. For the homologybased strategy putative protein sequences were searched against Drosophila melanogaster protein databases (FLYBASE. ORG) using the DEL-TABLAST algorithm, assigning the resulting GO term to the locus. For the predictive-based strategy the CATH assignments for large sequence datasets 127-129 was implemented.
PREQUAL V1. 02 132 [-pptype all] was used to filter nt sequences before the alignment, which was performed with MACSE V2. 03  98  and filtered subsequently with HMMCLEANER 133 and GBLOCKS V0. 91B 134 under a "relaxed" condition 135. A ML gene tree was then generated as implemented in IQ-TREE2 V2. 1. 3 COVID-EDITION 136 [--sampling GENESITE -m MFP --keep-ident] and used in the adaptive branch-site random effects likelihood (ABSREL) method  137, 138  as implemented in the HYPHY batch language  139. Enrichment of GOTERMS was performed using a combination of two different approaches, the HYPERGTEST algorithm, implemented in the GOSTATS package  140  for R [annotation org. Dm. eg. db; conditional TRUE; testdirection over), and GOATOOLS 141 (P value cutoff 0. 05); both using scOGs with a putative sign of diversifying positive selection (adjusted P value < 0.
The QS analysis was run using the scOG supermatrix described above, specifying quartet likelihood calculations with 100 replicates (i. e. , number of quartets draw per focal branch). Finally, Bayesian algorithm of MCMCTREE V4. 8A (from PAML package) 152 with approximate likelihood computation was used to estimate divergence times within the 64 Nymphalids. First, branch length by ML were estimated and then the gradient and Hessian matrix around these ML estimates were computed under MCMCTREE using the DNA supermatrix. As no fossils of Heliconiini or closely related tribes are known, we used four secondary calibration points (Supplementary Data 3; supplementary Fig. 21 ). These calibration points where selected based on their consistency with previous phylogenetic studies of Heliconiini and From the TIMETREE database  8, 153.
These tests were applied on all triplets extracted from scOG gene trees within Heliconiini, and the resulting P-values were then corrected for multiple testing using the Benjamini-Hochberg procedure with a false discovery rate (FDR) cutoff of 0. 05, using the provided R script BLT_DCT_FBRANCH_WRAPPER. R. We considered only introgression events that were significant in both tests. Because these tests are not independent, since different triplets may contain overlapping taxa, the correction results in more conservative tests 156 , and the inferred FDRs may be somewhat inaccurate due to the statistical correlation of some triplets (A. Suvorov personal communication). Dsuite 156 was then used to plot the results in a heatmap plot 14 using the script GETTRIPLETSSUPPORTINGINTROGEVENTS. PY available at the GITHUB page of Dsuite.
Furthermore, each OG was analyzed to contrast differential evolutionary pressures between Eueides and Heliconius species. Each OG of all these gene families was processed as follows: first, in the first step OGs subject to expansion (based on CAFE results) were excluded, alignments and gene trees were generated as in the signature of selection pipeline (see below). Second, from the tree, nodes with low bootstrap support (<0. 90) and only the putative in-paralog 160 with the shortest branch from its ancestor were retained (custom python script REMOVEINPARALOGFROM- TREE. PY) available at https://github. com/francicco/-ComparativeGenomicsOfHeliconiini. This procedure was implemented to generate near single-copy ortholog groups (nscOGs) to exclude paralog biases in follow-up analyses. NscOGs were then realigned separately using nt sequences for Eueides ssp.
and Heliconius ssp. , as in the signature of selection pipeline (see below). In each species group the mean ω was estimated (hyphy acd Universal $FASTA MG94CUSTOMCF3X4 Global) as implemented in  HYPHY. For each nscOG we tested if the gene experienced a differential selective regime testing for relaxation/intensification the Heliconius branches compared with the Eueides ones with RELAX 161 , as implemented in the HYPHY batch language 139 (P-value threshold 0. 05).
Annotating conserved non-exonic elements To identify conserved non-exonic elements (CNEEs), we used the 63way whole genome alignment and the PHAST V1. 4 package  162, 163. Using the E. isabella coding gene annotation we first extracted 4-fold degenerate sites using HAL4DEXTRACT with the flag -conserved; and extracted the alignment with HAL2MAFMP. PY with --noAncestors --noDupes flags, and used PHYLOFIT to estimate the initial neutral model. We reduced the dataset to a subset of 10 species representing the whole Heliconiinae phylogeny (S. mormoria, Dryadula phaetusa, Dione juno, E. isabella, H. erato, H. charitonia, H. sara, H. aoede, H. doris, H. melpomene), which significantly reduced computational time, making this step feasible. To assess the convergence of the model we performed five independent runs of PHYLOFIT using --init-random, the strand-symmetric version of the generalized time reversible model (SSREV), to check likelihood stability.