Linkage mapping and genome annotation give novel insights into gene family expansions and regional recombination rate variation in the painted lady (Vanessa cardui) butterfly

Genomics 114 (2022) 110481 Available online 14 September 2022 0888-7543/© 2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Linkage mapping and genome annotation give novel insights into gene family expansions and regional recombination rate variation in the painted lady (Vanessa cardui) butterfly Daria Shipilina a , b , * , 1 , Karin N¨ asvall a , 1 , Lars H¨ o¨ ok a , Roger Vila c , Gerard Talavera d , Niclas Backstr¨ om a a Evolutionary Biology Program, Department of Ecology and Genetics, Uppsala University, Norbyv¨ agen 18D, 75236 Uppsala, Sweden b Swedish Collegium for Advanced Study, Thunbergsv¨ agen 2, 75236 Uppsala, Sweden c The Butterfly Diversity and Evolution Lab, Institut de Biologia Evolutiva, Passeig Martim de la Barceloneta 37-49, 08003 Barcelona, Spain d Institut Bot` anic de Barcelona (IBB), CSIC-Ajuntament de Barcelona, Passeig del Migdia s/n, 08038 Barcelona, Spain ARTICLE INFO Keywords: Genomics Recombination Linkage map Gene family Painted lady Lepidoptera ABSTRACT Characterization of gene family expansions and crossing over is crucial for understanding how organisms adapt to the environment. Here, we develop a high-density linkage map and detailed genome annotation of the painted lady butterfly (Vanessa cardui) - a non-diapausing, highly polyphagous species famous for its long-distance migratory behavior and almost cosmopolitan distribution. Our results reveal a complex interplay between regional recombination rate variation, gene duplications and transposable element activity shaping the genome structure of the painted lady. We identify several lineage specific gene family expansions. Their functions are mainly associated with protein and fat metabolism, detoxification, and defense against infection - critical processes for the painted lady’s unique life-history. Furthermore, the detailed recombination maps allow us to characterize the regional recombination landscape, data that reveal a strong effect of chromosome size on the recombination rate, a limited impact of GC-biased gene conversion and a positive association between recombination and short interspersed elements. 1. Introduction The genomic era opens up opportunities for investigating relationships between genotypes and complex phenotypes on a novel level and for a better understanding of genome evolution. Combinations of different approaches can lead to novel insights into the dynamics of recurring duplications, deletions and other types of structural rearrangements, for example, by assessing molecular mechanisms and evolutionary consequences of gene family expansions and contractions, the activity of selfish genetic elements (e.g. transposable elements, TEs) and recombination rate variation. Gene duplication has since long been recognised as an important mechanism for generating novel genetic material for natural selection to act upon [41,79,120], and gene family expansions and contractions are important sources for generation of phenotypic diversity [24,55]. Comparative approaches, such as orthology analysis, allow for identification of expanding or contracting gene families and annotation of orthogroups with functional relevance in the evolution of lineagespecific traits. This approach might be beneficial for investigating complex phenotypes such as migratory behavior or polyphagy, where combined effects of different types of genetic changes likely underlie the trait [92]. Since the spearheading work by McClintock [72], transposable elements (TEs) have been acknowledged as major contributors to different types of evolutionary change in eukaryotes [52,54]. Transposable elements are capable of self-replication within the host genome and can mediate both small scale deletions and duplications, large scale chromosomal rearrangements [54] and considerable genome size expansions [82,102]. In addition, TE insertions can affect gene function when regulatory or coding regions are targeted. Therefore, characterization of the TE repertoire is key to understanding the microevolutionary dynamics within the genome of a species and the potential effects of TE activity on * Corresponding author at: Evolutionary Biology Program, Department of Ecology and Genetics, Uppsala University, Norbyv¨ agen 18D, 75236 Uppsala, Sweden. E-mail address: [email protected] (D. Shipilina). 1 Authors contributed equally. Contents lists available at ScienceDirect Genomics journal homepage: www.elsevier.com/locate/ygeno https://doi.org/10.1016/j.ygeno.2022.110481 Received 4 May 2022; Received in revised form 1 September 2022; Accepted 10 September 2022 Genomics 114 (2022) 110481 2 trait variation within populations and between species. Besides gene duplication and TE activity, recombination is a process crucial for evolutionary innovation. Meiotic recombination shuffles existing segregating genetic variants, resulting in the generation of novel haplotypes [81]. Recombination also influences selection efficiency by directly preventing the accumulation of deleterious alleles (Müller’s ratchet) and breaking the physical linkage between mutations with different selective effects (Hill-Robertson effects). The rate of recombination can vary on different scales. Of particular interest for population genetic processes is the variation in recombination rate between different genomic regions. Such spatial variation in the recombination rate has been observed in many different organisms [98,103]. However, besides detailed recombination maps in the butterfly genus Heliconius [71], little is known about how the rate of recombination rate varies across chromosome regions in Lepidoptera and how recombination is associated with different genomic features [37,101]. As indicated above, incorporating different approaches is essential for studying the genetic underpinnings of complex phenotypes and the mechanisms governing microevolutionary processes. The painted lady, Vanessa cardui, represents a key study system for a wide array of evolutionary studies. It is the most wide-spread of all butterfly species [99], and its migratory behavior includes a diverse repertoire of distinct phenotypes. In general, migratory butterfly species need to sustain longdistance flight and have well developed navigational abilities [23,35]. Therefore, traits related to energy metabolism, sensory reception and the flight machinery have likely been under strong directional selection. In contrast to many other migratory butterflies inhabiting temperate zones, the painted lady is a non-diapausing, multigenerational migrant, with an annual migratory circuit covering areas with extreme environmental heterogeneity [73,100]. Despite the high risks associated with such a migratory lifestyle, painted ladies have successfully colonized almost all continents, and the species harbors high levels of genetic diversity, indicating a large effective population size [33]. This could also be a consequence of the species’ ability to utilize a wide range of host plants [2,21]. Until the era of high-throughput sequencing, the possibilities to gain insights into how the migratory and generalist lifestyle has been manifested at the level of the genome have been limited: genetic basis of migratory behavior in insects has only been investigated in a few model species (e.g. the monarch butterfly, Danaus plexippus) so far [74]. A key step for genomic analyses is the development of a highcontiguity genome assembly of the focal species and a thorough genome annotation. A powerful method to ensure the spatial correctness of a chromosome level physical assembly is construction of a linkage map. In this study, we present the first detailed linkage map of the painted lady and verify scaffolds from a previously available genome assembly based on long-read sequencing technology [67]. We use the genome annotation and linkage information to quantify lineage-specific patterns of gene family evolution, relative TE abundance and how the regional recombination rate variation is associated with genomic features in the painted lady. Our analyses complement earlier efforts to establish genomic tools for this species [26,121] and give novel insights into the overall genome structure, recombination rate variation and lineage-specific gene family expansions in this species, information that informs on the molecular mechanisms underlying genome evolution in butterflies in general and the formation of the complex migratory phenotype and generalist lifestyle of the painted lady in particular. 2. Results 2.1. Linkage map and genome annotation To verify a chromosome level assembly of the painted lady [67] and to get access to detailed recombination rate data, we constructed a pedigree-based linkage map. The total distance of the linkage map was 1516 centiMorgan (cM) and contained 1323 markers. When anchored on the 424 Mb physical assembly, the average marker density was 3.09 markers / Mb. The genome assembly was highly collinear with the marker order in the linkage map (Pearson’s correlation coefficient; R = 0.91–1.00, p-value >1.00 ×10 −04 , Fig. S1) and consisted of 30 autosomes and the sex chromosomes Z and W. The high collinearity between linkage groups and assembled scaffolds, the large scaffold N50 (14.6 Mb) and high BUSCO scores (97% complete arthropod genes) confirm that the scaffolds in the assembly essentially represent complete chromosomes that could be used for accurate characterization of genomic features and quantification of regional recombination rate estimates. In total, TEs constituted >150 Mb (37.40%) of the assembly and LINEs and SINE were the most abundant of the characterized repeat classes (Table 1). After automatic annotation and subsequent manual curation, 13,161 protein-coding genes were identified (including 89.90% BUSCO genes), of which 12,209 had functional annotation information (Table 1). Visual inspection of the spatial distribution of genes and TEs along chromosomes revealed rather similar distributions of repeat classes between autosomes and the Z-chromosome, but also an observable excess of repeats on smaller autosomes and a striking difference in repeat composition and gene density on the W-chromosome (Fig. 1). 2.2. Synteny The level of large-scale structural conservation of the painted lady genome was assessed by comparing gene order on the painted lady chromosomes to two previously available high-contiguity lepidopteran genome assemblies positioned at different levels of divergence in the lepidopteran tree of life, the silkmoth (Bombyx mori) and the postman butterfly (Heliconius melpomene). Overall, the synteny was highly conserved between the painted lady and the other species, but chromosomes 28 and 26 mapped to the same chromosome (24) in B. mori and the previously described fusions of several chromosomes in the H. melpomene genome [27] could also be verified (Fig. 2). In summary, this confirms that the painted lady karyotype is highly similar to the inferred ancestral butterfly karyotype [3]. 2.3. Gene family evolution To investigate the turnover of specific gene families in the painted lady, we analyzed a set of nine representative nymphalid species with detailed annotation information (see methods). The non-migratory Kamehameha butterfly (Vanessa tameamea) was included to assess Table 1 Linkage map, genome assembly and annotation statistics. Linkage map Total map length (cM) 1516 Number of markers 1323 Markers per physical distance (N / Mb) 3.09 Genome assembly Scaffold N50 (bp) 14,615,999 GC content 33.41% Total repeat proportion 37.40% Repeat content (% of total repeat proportion) SINEs 7.30% LINEs 14.94% LTR elements 2.47% DNA elements 3.04% Simple and unknown repeats 30.17% Gene annotation BUSCO genes 89.90% Number of protein coding genes 13,161 Number of genes with functional annotation 12,209 D. Shipilina et al. Genomics 114 (2022) 110481 3 Fig. 1. Distribution of repeat classes and genes as estimated along the painted lady chromosomes (100 kb windows). Density (% of the window covered) of different TE classes are illustrated with distinct colors cumulatively added on top of each other above the X-axis and density of genes below the X-axis (legend to the top right). D. Shipilina et al. Genomics 114 (2022) 110481 4 differences in gene family evolution between sedentary and migratory lineages within the Vanessa genus. We found that 93.2% (1,288,332) of the total number of genes from the nine nymphalid species were clustered in 14,027 orthogroups. The percentage of genes assigned to orthogroups varied from 86.7 to 99.6% in the different species (Table S1). In the painted lady, 96.4% (12,692) of the annotated genes were assigned to 10,361 orthogroups with 19 lineage-specific orthogroups containing 63 genes (Table S1). Within the Vanessa genus, 65 expansions had occurred on the ancestral branch, 648 on the V. cardui branch and 1563 on the V. tameamea branch. We used a maximum likelihood model to detect genes with distinct gene family expansion rates in the painted lady compared to the other species. The analysis showed that 12 orthogroups were significantly expanded in the painted lady. These orthogroups contained 77 genes, of which 34 had associated GO-terms (Fig. 3). Among the largest expanded gene families were two classes of proteases, a lipoprotein receptor and the Lepidoptera-specific moricin immune-gene family. Analysis of the spatial distribution of extended orthogroups revealed clustering/tandem duplications for all except one of the orthogroups (Fig. 3C). Significantly enriched GO-terms for expanded gene families in the painted lady were predominantly associated with protein degradation, muscle function and development, and fatty acid and energy metabolism (Fig. 3A). Multiple ontology terms were shared between expanded orthogroups, pointing towards similar functions associated with the different gene families (Fig. 3B). Additionally, we identified gene families with a distinct gene expansion rate in both the painted lady and the monarch butterfly Danaus plexippus - the latter a key model organism for insect migration studies - and compared those to the other nymphalids. This analysis revealed 11 orthogroups with a higher expansion rate and 29 orthogroups with genes specific to these two lineages. The common orthogroups included 112 genes and were significantly enriched for GOterms predominantly associated with metabolic processes, defense against infection and neuronal activity (Fig. S2). 2.4. Patterns of recombination rate variation 2.4.1. Global and chromosome specific recombination rates The development of a detailed linkage map allowed both for estimating the global recombination rate in the painted lady and to investigate potential regional recombination rate variation and association with genomic features. The average, genome-wide recombination rate was 3.81 cM / Mb (W-chromosome excluded), but there was considerable inter-chromosomal variation (2.21–8.00 cM / Mb; Table S2, Fig. S3), with a significantly higher rate on shorter chromosomes than on longer chromosomes (Spearman’s rank correlation, ρ = − 0.83, pvalue =6.51 ×10 −07 ;Fig. 4). The recombination rate on the Z-chromosome was 3.09 cM / Mb, lower than the average unweighted autosomal rate. However, the recombination rate on the Z-chromosome was not lower than expected given the overall negative correlation between recombination rate and chromosome size (Fig. 4A). Besides the negative association between chromosome size and recombination rate, we also found significant negative associations between chromosome size and GC-content ( ρ = − 0.65, p-value =8.35 ×10 −05 ) and repeat density ( ρ = 0.77, p-value =1.37 ×10 −06 ), and a positive association with gene density ( ρ =0.68, p-value =2.63 ×10 −05 )(Fig. 4B-D). 2.4.2. Intra-chromosomal variation in recombination rate and associations with genomic elements To quantify potential regional variation in recombination rate within chromosomes, we estimated the recombination rate in 2 Mb nonoverlapping windows along each individual chromosome. The average rate across windows was similar to both the global rate estimate across chromosomes (4.05 +/−2.45 cM / Mb) and the overall chromosome level estimates (2.58–7.53 cM / Mb, W-chromosome excluded). The recombination rate estimates for individual windows ranged between 0 and 14.79 cM / Mb (Fig. S3, Table S2) and visual inspection revealed a bi-modal distribution with reduced recombination rate in the center of chromosomes and towards chromosome ends (Fig. 4 E-H). To test this observation formally, we analyzed the difference in recombination rate between bins representing five relative distance intervals from the center of the chromosome for all chromosomes combined and found that Fig. 2. Synteny between the painted lady and a) the silkmoth (Bombyx mori) and b) the postman butterfly (Heliconius melpomene) chromosomes, respectively. The painted lady autosomes are sorted and named according to length and the sex chromosomes are indicated with Z and W. Orthologous genes are connected with lines. Colors represent individual chromosomes. D. Shipilina et al. Genomics 114 (2022) 110481 5 the recombination rate was significantly lower in the center (first bin), significantly higher in the flanking terminal (fourth) regions and then again lower at the terminal end (Wilcoxon rank sum tests, p-values = 3.70 ×10 −02 - 6.30 ×10 −15 ; Fig. S4). To assess potential relationships between the recombination rate and genomic features in more detail, we first investigated different associations between the window-based recombination rate estimates and variation in nucleotide composition and proportions of different TEs and genes. The W-chromosome was excluded from this analysis since it is non-recombining in Lepidoptera. We found that the GC-content increased towards the ends of chromosomes and was positively associated to the regional recombination rate ( ρ =0.32, p-value =3.68 × 10 −06 ). Gene density was homogeneous across chromosomes, with only a minor increase towards the chromosome center, and was negatively associated with the recombination rate ( ρ = − 0.19, p-value =7.27 × 10 −03 ). We found a significant positive association between the overall repeat proportion and the recombination rate ( ρ =0.35, p-value =3.48 ×10 −07 ; Fig. 5), and this pattern was consistent for all repeat classes, but strongest for SINEs ( ρ =0.42, p-value =2.63 ×10 −10 ) and weakest for LTRs ( ρ =0.14, p-value =4.04 ×10 −02 ). The association between recombination rate and proportion of LTRs was, however, not significant when only including autosomes ( ρ =0.11, p-value =11.04 ×10 −02 ; Fig. 5, Fig. S5). To disentangle the relative strength of associations between the regional recombination rate and genomic features, a multiple linear model was implemented with recombination rate as the dependent variable. As explanatory variables we used chromosome length, chromosome type, GC-content, proportion of genes (CDS) and proportions of all different classes of TEs. We found that the regression model was significant (df =197, F =7.73, p-value =5.36 ×10 −09 ) and explanatory variables in the model accounted for 21% of the variation in recombination rate (R2 =0.24, adjR2 =0.21). Most of the variation was explained by the positive association with the proportion of SINEs (Estimate 1.59, p-value =9.75 ×10 −04 ) and the negative association with chromosome size (Estimate −0.48, p-value =4.88 ×10 −02 ; Fig. 5, Table S3). Fig. 3. A) Significantly (p-value <0.05 after FDR-correction) enriched gene ontology (GO) terms associated with expanded gene families in the painted lady. The bars show the number of genes associated with each GO-term. The different GO-categories are biological process (BP), cellular compartment (CC) and molecular function (MF). B) Orthogroups with significantly enriched GO-terms. Shared GO terms between ontology terms (biological process category only) are shown in connecting lines. C) Spatial distribution of genes from extended orthogroups identified in BadiRate analysis. D. Shipilina et al. Genomics 114 (2022) 110481 6 Finally, we explored whether gene expansions could be associated with other genomic features, and we therefore compared TE abundance in the regions with and without gene gains. The mean densities of LTRs, LINEs and DNA transposons were higher in regions with gene gains (Wilcoxon rank sum test, p-value 3.1 ×10 −03 - 6.0 ×10 −04 ; Fig. S6), as was mean GC-content (p-value 3.0 ×10 −02 ). The gene densities or recombination rates did not differ between regions with or without gene gains (Wilcoxon rank sum tests, p-values =9.1 ×10 −01 - 8.0 ×10 −01 ; Fig. S6). 3. Discussion 3.1. The genome of the painted lady butterfly Here we present detailed results on the genomic architecture and regional recombination rate variation in the painted lady. The data paves the way for understanding the interplay between molecular mechanisms and micro-evolutionary processes shaping the genome of butterflies in general and provide the first insights into the links between genomic features and the unique lifestyle of this species. The rapid technological advances and dropping costs of DNA-sequencing methods have led to a staggering development rate of high-quality genome assemblies, including many butterfly species [20,34,63,96,116], and the availability of genomic resources will probably increase almost exponentially in the near future, as a result of the Darwin tree of Life (http s://www.darwintreeoflife.org/), the European Reference Genome Atlas (ERGA; https://www.erga-biodiversity.eu/) and other similar initiatives. However, detailed and curated genome annotation data are more time-consuming and expensive to generate and therefore still limiting comparative/population genomic and genotype-phenotype association approaches, not the least in butterflies [27,42,106]. Another limiting factor for understanding both genome architecture in general, the relative effects of random and selective forces on sequence evolution and maintenance/loss of genetic diversity is that detailed recombination rate data are both laborious and time-intensive to gain, especially for natural populations. As a consequence, high-density recombination maps are still lacking for the vast majority of wild species where genome assemblies are now available. The detailed annotation information and the high-density linkage map for the painted lady developed here, therefore provide opportunities for both comparative studies on genome structure organization, population genomicand micro-evolutionary investigations in the entire Lepidoptera clade. Chromosome numbers have been shown to vary considerably between different butterfly and moth species; the haploid chromosome counts range from 5 to 223 [28,68]. In agreement with previous data [121], both the linkage map and the DToL genome assembly clearly showed that the painted lady has a total haploid chromosome count of 31. We confirmed high levels of synteny and gene order collinearity between the painted lady and the silkmoth, and the lineage specific chromosome fusions characterized before in the postman butterfly [27]. Hence, similar to other nymphalid butterflies, the painted lady has retained the inferred ancestral lepidopteran karyotype [3]. The annotation procedure revealed that the painted lady harbors a gene set (n = 13,161) close to the suggested core set in Lepidoptera [22,60] and a relatively low overall TE content. However, the TE content was significantly higher and the gene density lower on smaller chromosomes. 3.2. Sex chromosomes The assembly and annotation of sex-chromosomes, especially the non-recombining parts of sex-limited chromosomes (i.e. the W-chromosome in Lepidoptera), can be technically challenging due to the high density of repetitive elements. Up to date, there are only a few Lepidoptera species where the W-chromosome has been assembled and annotated [75]. Given the high-quality assembly we had access to, we performed annotation and manual curation of TEs and coding genes for Fig. 4. Top row: Associations between chromosome length and A) recombination rate, B) base composition, C) repeat and D) gene proportions. Chromosome length is given in megabases (Mb). Bottom row: Regional distribution of the recombination rate (E), base composition (F), repeat (G) and gene (H) density in 2 Mb windows along the chromosomes. All chromosomes were analyzed jointly and the x-axis shows the relative position (proportion of chromosomal length) from the center of the chromosomes. D. Shipilina et al. Genomics 114 (2022) 110481 7 the painted lady W-chromosome. In contrast to previous annotation [67], we could not confirm the presence of any protein coding genes. Gene models created on the preliminary annotation were not confirmed after manual curation and functional domain annotation. A lack of protein coding genes on the W-chromosome has also been observed in the silkmoth [1,75]. This apparent complete loss of protein coding genes on the Lepidoptera W-chromosome is obviously a consequence of the degradation process that has been well described for non-recombining parts of sex-chromosomes in many systems [10]. While having a size equal to an average autosome, the W-chromosome also demonstrated a significantly higher overall proportion of TEs, a larger fraction of longer TEs, and a different distribution of repeat classes compared to other chromosomes. Similar to the silkmoth and julia heliconian (Dryas iulia), the Wchromosome in the painted lady had a significantly higher proportion of LTRs and LINEs [59,75]. The proportion of SINEs was however much smaller on the W-chromosome than on the autosomes and the Z-chromosome. The higher accumulation of TEs is also an expected consequence of recombination suppression and comparatively low effective population size (N e ) of the W-chromosome (1/4 of the autosomes at equal sex-ratios), both as a consequence of Müllers ratchet and since the overall efficiency of selection against TE insertion is reduced for nonrecombining chromosomes [10]. The Z-chromosome is generally highly conserved in Lepidoptera [32] and it is the largest of all the painted Fig. 5. Correlation between recombination rate and density of genomic features. A) Summary of the linear model with regional recombination rate as the response variable. Each explanatory variable in the model is listed along the Y-axis and the relative estimated effect (X-axis), and error intervals are indicated with horizontal bars for each variable. B - H) Associations between the regional recombination rate and specific genomic features. Linear regression lines, Spearman’s correlation coefficients ( ρ ) and the corresponding ρ -values for significant analyzes are given. Gray dots/lines indicate autosomal regions and turquoise dots/lines Z-chromosome linked regions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) D. Shipilina et al. Genomics 114 (2022) 110481 8 lady’s chromosomes. We did not find any significant differences in gene or TE content on the Z-chromosome compared to the autosomes. 3.3. Gene family analysis Gene family expansions can provide the raw material for both neoand sub-functionalizing evolutionary directions, and the rate of gene duplication can be significantly higher than the rate of function-altering single nucleotide mutations [65]. However, most gene duplication events are probably deleterious [66] or effectively neutral, leading to a low probability of fixation of novel gene copies [30]. We found a comparatively low proportion of lineage-specific gene duplications in the painted lady, which could be a consequence of the large N e of the species [33], which translates to efficient selection against slightly deleterious variants. The majority of the significant gene expansions in the painted lady lineage clustered on single chromosomes - only a single gene family had expanded and dispersed across multiple chromosomes - suggesting that unequal crossing over has been the main mechanism behind gene family expansions. The painted lady has an extraordinary life-history and has become a quickly uprising complementary model organism for studying insect migration. Over most of the almost cosmopolitan distribution range [93], the painted ladies complete a multigenerational migratory circuit, where single individuals can migrate >4,000 kilometers during lifetime (Talavera and Vila, n.d.). In contrast to other migratory butterflies like the monarch and the red admiral (Vanessa atalanta), the painted lady is non-diapausing [93]. The genetic underpinnings of migratory behavior have only been preliminarily characterized for a handful of insect species [49,122] and have not been studied in painted lady before. The dissection of potential associations between genetic (and epigenetic) variants and complex phenotypes like migratory behavior requires a combination of multiple approaches. As the first step to understanding lineage-specific characteristics of the painted lady, we here focused on gene family evolution. Our results showed a limited number of genes with significant copy number expansions unique to the painted lady lineage. The expanded gene families were mainly associated with functions related to the transport of fatty acids, protein metabolism, and muscle structure and activity. Since migratory insects mainly use fat as an energy resource during migration [58,76,97,110], both the capacity to build up fat deposits and efficient sequestration of fatty acids have likely been under strong selection in the painted lady. Likewise, enhanced muscle structure and function should be advantageous for long-distance migrants compared to sedentary species. Therefore, efficient fine-tuning and optimization of fatty acid metabolism and increased muscle sustainability during migration could have been aided by the expansion of specific gene sets involved in those processes. Long-range migrants benefit from utilizing a multitude of different host plants since they will encounter dramatically different habitats, both during the lifespan of single migratory individuals and between consecutive generations. In contrast to the monophagous monarch butterfly, the painted lady can utilize >300 different larval host-plants in 11 plant families [2,21,78]. Two of the significantly expanded gene families in the painted lady (UDP-glycosyltransferase, carboxylesterase) were associated with polyphagy and detoxification [18,39,77]. The UDP-glycosyltransferase superfamily includes Lepidoptera-specific subfamilies associated with a variety of functions, such as affinity for plant secondary metabolites [44,69]. In the painted lady larvae, one UDP-subfamily is upregulated in response to utilization of an extended range of hostplants [21]. Copy-number expansions of these detoxifying gene families could have allowed the painted lady to increase the range of host plants that can be utilized and consequently paved the way for developing the non-diapausing, multigenerational, long-distance migratory lifestyle. The wide range of habitats that long-distance migratory species encounter also probably means that they are exposed to many more different pathogens than sedentary species. Our analysis revealed that the Lepidoptera-specific gene moricin, associated with inducible antimicrobial peptides [38], was significantly expanded in the painted lady. An increase in the number of moricin copies could have increased the efficiency of defense against a larger suite of pathogens. Previous investigation of the genetic basis of migratory behavior in the monarch butterfly identified candidate genes associated with orientation, chemoreception and regulation of the circadian clock [119,122]. Migratory behavior has evolved independently multiple times within the Papilionoidea clade [25] and in the Vanessa genus [108], and the life histories of the monarch butterfly and the painted lady are distinct. However, long-distance migration should put selective pressure on similar traits (e.g. navigation, energy metabolism, muscle endurance), and it is therefore possible that specific gene categories have been under selection in independent lineages. Significantly expanded gene families shared between the painted lady and the monarch were enriched for functions associated with various metabolic processes, defense against pathogens and neuronal activity, all of which can be associated with migratory behavior. One gene family with an especially pronounced expansion was vacuolar ATPases, ATP-dependent proton pumps involved in membrane ion transport [113]. Given the unique expansion of this gene family in both species, we speculate that copy number increase could be involved in flight muscle coordination and/or ion transport for maintenance of homeostasis during long periods of flight. In this study, we get a first glimpse of the specific genes that have undergone copy number expansions in the painted lady specifically and independently in the two migratory species. The functions associated with the expanded gene families can be coupled to the evolution of longdistance migratory behavior. However, further studies of independent migratory and sedentary sister species, in combination with detailed population genetic analysis and functional verification will be necessary to dissect the genetic underpinnings of migratory behavior in butterflies in detail. 3.4. Patterns of recombination rate variation Detailed data on recombination rate variation are crucial for understanding the relative effects of genetic drift and selection on levels of genetic diversity. Understanding how recombination breaks down linkage disequilibrium is also important for association studies aimed at coupling genetic variation to phenotypic traits. Despite their importance, detailed recombination maps are only available for a handful of butterfly species [14,20,27,89,96,105]. In some butterfly species, linkage maps have been used to improve and/or verify the correctness of physical genome assemblies, but the recombination rate has not been assessed. Here we developed a high-density linkage map based on segregation information in a pedigree with 95 offspring. The map contained >1,300 ordered markers and the overall density was >3 markers per Mb. Despite being based on a single pedigree, the genetic map developed here revealed a recombination landscape in strong agreement with what has been observed in other butterflies [27,71]. This indicates that the painted lady genetic map accurately reflects the historical recombination landscape in the species. We estimated the genome-wide average recombination rate in the painted lady to be 3.81–4.05 cM / Mb, dependent on the method applied. The global rate was in the lower end of recombination rate estimates from other Lepidoptera species, which have been in the range from 2.97 to 4.0 cM / Mb in the silkmoth [115,117] to 5.5–6.0 cM / Mb in different Heliconius species [47,104]. We found a significant negative association between chromosome length and the recombination rate in the painted lady. This is a consistent pattern found across many organism groups and likely a consequence of that at least one crossover event is necessary for correct segregation of chromosomes during meiotic division in the recombining sex, leading to a higher recombination rate per unit length for shorter chromosomes [37,50,71]. D. Shipilina et al. Genomics 114 (2022) 110481 9 Butterflies and moths have holocentric chromosomes, i.e. they lack distinct centromere regions, which might lead to an expectation of a uniform distribution of recombination events. In the painted lady we observed a bimodal distribution of recombination events along chromosomes, with an increased recombination rate away from the center and significant drops at the chromosome ends. This distribution is in agreement with previous observations, both in Lepidoptera and in other animals with different centromere types [37,71]. A possible explanation for this pattern is mechanical or tension interference between chiasmata when >1 recombination event occurs on the same chromosome [37]. However, in the holocentric Caenorhabditis elegans, the number of recombination events is limited to precisely one per chromosome per meiosis, but there is still a strong bimodal pattern of recombination rate variation along chromosomes in this species [12]. An alternative explanation could be that synaptonemal complexes are directed towards the flanking regions, when the telomeres attach to the nuclear wall [91]. The reduced recombination rate at chromosome ends is also consistent with earlier observations and could potentially be attributed to selection against synaptonemal complex formation at chromosome ends, due to a higher risk of ectopic recombination in these generally repeat-rich regions [95]. Since recombination is directly associated with the efficacy of selection, a negative correlation between the regional recombination rate and number of repeats would be expected if TE insertions predominantly are deleterious. Such associations have been observed in many organisms, although the relationship between TE-abundance and the recombination rate varies to some extent across species and different TEclasses [53,88]. In the painted lady, we observed a significant positive association between TE-abundance and the regional recombination rate, predominantly driven by a strong effect of SINE density. An explanation for the strong association between SINE density and recombination rate could be SINE-mediated recombination, as has for example been described in humans [29]. In addition to the strong positive association between SINE density and recombination rate on the autosomes and the Z-chromosome, the absence of SINEs on the non-recombining W-chromosome supports that SINEs might be able to hijack the recombination machinery. However, we can not exclude other factors affecting both the recombination rate and the proliferation efficiency of SINEs. For example, both synaptonemal complexes and SINE insertions might be directed towards regions of more open chromatin structure. In contrast with results from similar studies in other organism groups [9,51], we observed a negative association between the recombination rate and gene density. This is likely a consequence of the strong association between recombination rate and chromosome size, since the association with gene density was insignificant when chromosome size was included as an explanatory variable. The observed weak positive association between GC-content and recombination is in agreement with the limited effect of GC-biased gene conversion (gBGC) in butterflies [17]. We did not find any association between recombination rate and the presence of extended orthogroups, which would be expected if gene duplication is associated with unequal crossing-over. This could possibly be a consequence of the more efficient removal of deleterious duplications in regions with higher recombination rate. However, repetitive elements can trigger ectopic recombination which can explain the observed significant positive association between gene gains and density of LTRs, LINEs and DNA elements in the painted lady. 4. Conclusions In this study, we present detailed annotation and recombination rate information for the painted lady butterfly (Vanessa cardui), a species with remarkable life-history traits such as long distance migration, continuous direct development and a capacity to utilize many different types of larval host plants. We analyzed lineage-specific gene family expansions and found that expanded genes were mainly associated with fat and protein metabolism, detoxification and defense against pathogens. A detailed TE-annotation revealed that several TE-classes were positively associated with the presence of gained genes, potentially indicating their involvement in ectopic recombination. Recombination rate variation was negatively associated with chromosome size and positively associated with the proportion of short interspersed elements (SINEs). We conclude that the genome structure of the painted lady has been shaped by a complex interplay between recombination, gene duplications and repeat activity and provide the first set of candidate genes potentially involved in the evolution of migratory behavior in this almost cosmopolitan butterfly species. 5. Methods 5.1. Linkage map 5.1.1. Sampling and DNA-extraction Offspring from one painted lady female were reared on thistles (Cirsium vulgare) in the greenhouse until pupation. The bursa copulatrix of a female was examined and only one spermatophore was detected, indicating that a single male had sired all offspring. The offspring were snap frozen in liquid nitrogen and stored in −20 ◦C until DNA extraction. DNA was extracted from thorax tissue of the female and an abdominal segment of the offspring pupae, using a modified high salt extraction method [7]. The quality of the DNA was analyzed with Nanodrop (ThermoFischer Scientific) and the yield was quantified with Qubit (ThermoFischer Scientific). Extracted DNA was digested with the restriction enzyme EcoR1 according to the manufacturer’s protocol, using 16 h digestion time (ThermoFischer Scientific). DNA fragmentation was verified with standard gel electrophoresis. Digested DNA from 95 offspring with the highest yield and the dam was shipped to the National Genomics Infrastructure (NGI, see acknowledgements) in Stockholm for library preparation (standard protocol), individual barcoding and multiplex sequencing using 2 ×151 bp paired-end reads on one NovaSeq6000 S4 lane. 5.1.2. Building the linkage map The quality of the raw reads was assessed with FastQC [8]. The reads were filtered using the Stacks2 modules clone_filter to remove PCRduplicates and process_radtags to filter for quality. We evaluated phredscore in sliding windows covering 15% of the read length and removed reads with mean score below 10 [19]. Removal of reads with unassigned bases and truncation to 125 bp was done using option -c, and –disable_rad_chec was applied to keep reads with incomplete RAD-tags. We mapped the filtered reads to the previously published genome assembly [67] using the bwa mem algorithm [61] with default options. Resulting bam files were sorted with samtools sort [62] and filtered with samtools view −q 10 (only reads with mapping quality score above 10 were retained). A custom script was applied to retain reads with unique hits only. The mapping coverage was analyzed with Qualimap [80]. The offspring were defined as females if the coverage on the Zchromosome was <75% of the average coverage over all chromosomes and as males if the coverage was >75%. Samtools mpileup was used for variant calling using minimum mapping quality (−q) 10 and minimum base quality (−Q) 10 [62]. The variants were then converted to likelihoods with Pileup2Likelihoods in LepMap3 using default settings [85]. The LepMap3 protocol [85] with some modifications was used to construct the linkage map (Supplementary methods 1). 5.2. Genome annotation and whole genome statistics 5.2.1. Genome assembly statistics With very few exceptions, the order of markers in the linkage map was in agreement with the physical order in the assembly. We therefore did not make any corrections to the physical assembly before further analysis. Standard genome assembly summary statistics were calculated D. Shipilina et al.