Demographic history inferences

The demographic trajectories, inferred separately for each chromosome, jointly revealed the demographic history of the Swedish L. sinapis population. These trajectories shared three main features: firstly, a maximum effective population size (N_e) of ≈10⁶ approximately 10,000 generations before present (BP), preceded by a slight increase matching with the end of the Last Glacial Period (LGP). More recently, N_e gradually declined until stabilizing about 100 generations BP. Contemporary estimates of N_e oscillate between 10³ and 2.5 × 10⁴ for the different chromosomes (Fig. 1).

Figure 1.

Demographic trajectories inferred for the different chromosomes in the Swedish population of the wood white. Six out of the ten individuals with the highest sequencing coverage were used to represent the entire geographic distribution. Only the 100,000 most-recent generations of the trajectories are shown, as all useful variation in N_e happens within that range. The number of generations is equivalent to the identical number of years only if we assume univoltinity over time. The black line indicates the mean N_e across chromosomes.

Open in new tabDownload PowerPointLink to figure

To test for consistency in the demographic inference, the analyses were repeated using another six randomly selected distinguished individuals. There was no significant difference in N_e (mean difference = 1618 across chromosomes) between the two analyses (Wilcoxon's test, V = 242, P = 0.38), indicating that the selection of distinguished individuals has a minor effect on the inferred demographic trajectories.

Recombination rate variation and distribution

The genome-wide recombination rate was 7.37 cM/Mb, with individual chromosomes ranging between 3.5 and 15.3 cM/Mb. There was a trend toward lower recombination rates on larger chromosomes, but it was not significant (Spearman's rho = −0.292, P = 0.06; Supplemental Fig. 1). Also, the recombination rates on autosomes (7.65 cM/Mb) and Z Chromosomes (7.03 cM/Mb) were not significantly different (Wilcoxon's test, W = 41, P = 0.92).

Visual inspection of the recombination maps revealed a considerably reduced rate toward chromosome ends (Fig. 2A). A formal analysis showed that the recombination rate in the subtelomeric regions (five 100-kb windows at each chromosome end, respectively) was significantly reduced (2.46 cM/Mb) compared to proximal chromosome regions, that is, outside of the subtelomeric regions (Wilcoxon's test, W = 1,310,856, P-value = 3.18 × 10⁻⁶⁶). This, however, does not necessarily reflect a low recombination rate in the telomeres per se, because these were not covered in the wood white genome assembly (Höök et al. 2023). We found a significant reduction in recombination rate in 50 out of the 58 subtelomeric regions (29 chromosome pairs, 2n = 58; Supplemental Table 1). Four out of the eight chromosome ends with no significant reduction in recombination rate match observations of segregating fusion/fission polymorphisms; these involve chromosome pairs 18 and 25, 11 and 26, and 5 and 27 (Höök et al. 2023).

Figure 2.

Recombination rate estimates and hotspot characterization. (A) Estimates of the recombination rate (cM/Mb) in 1-Mb windows, chromosomes ordered by chromosome length. (B) As an example, we show a detailed recombination map for Chromosome 24. The inferred regional recombination rate is represented by gray bars and the orange line represents the threshold for calling recombination hotspots, equivalent to 10× the background recombination rate for a given interval. Red lines underneath the plot indicate the positions of inferred recombination hotspots.

Open in new tabDownload PowerPointLink to figure

Besides the reduced recombination rate in subtelomeric regions, we found other differences in regional variation between chromosomes. First, some chromosomes showed a unimodal recombination rate distribution, with higher values in central regions and a progressive decrease toward the ends. This was a frequent (but not exclusive) pattern in shorter chromosomes (<20 Mb, Fig. 2A). In contrast, the recombination landscape was bimodal along most longer chromosomes (>20 Mb, Fig. 2A), with a central region of reduced recombination in addition to the reduction in subtelomeric regions.

Recombination hotspot and coldspot identification

Informed by coalescent simulations, we developed thresholds for the identification of recombination hotspots and coldspots, taking into account the population's demographic history and the potential stochasticity of our method of recombination rate inference (see Methods for further details). Based on these thresholds, 3124 recombination hotspots were classified (Fig. 2B). These were on average 1656 bp long, with a mean recombination rate of 94.1 ± 62.5 cM/Mb (i.e., an ∼13-fold increase over the genome-wide estimate) and a maximum rate of 708 cM/Mb. This translates to 10.7% of historical recombination events in the population occurring within recombination hotspots, despite them constituting only 0.87% (5.2 Mb) of the genome. Recombination hotspots were found at a significantly lower frequency in the terminal 10% regions of the chromosomes (5% on each end; Wilcoxon's W = 220.5; P = 4.94 × 10⁻⁴; Supplemental Fig. 2A). Permutation analysis revealed significantly lower long terminal repeats (LTRs) (P < 0.001, Fig. 3B) and LINE densities (P < 0.001, Fig. 3D), and a higher short-interspersed nucleotide element (SINE) density (P < 0.001, Fig. 3C) in the hotspots, although DNA transposon density did not deviate significantly (P = 0.51, Fig. 3A) from the genome wide average. No functional gene categories or specific sequence motifs were significantly enriched in hotspots (P > 0.05, Supplemental Table 2).

Figure 3.

Occurrence of transposable elements (TEs) in hotspots. The histograms show the density distributions of (A) DNA, (B) LTR, (C) SINE, (D) LINE, for 3124 randomly sampled 1656-bp genomic windows (number of, and average length of identified hotspots). The red lines indicate the observed mean TE densities in the defined hotspot regions.

Open in new tabDownload PowerPointLink to figure

We also identified 1283 recombination coldspots, that is, regions with considerably reduced recombination rate as compared to the genomic average (see Methods). The average length of the coldspots was 30 kb, with 70 of them being longer than 100 kb. Despite the lower frequency of coldspots than hotspots, they represented a larger proportion of the genome (38.2 Mb, 6.37%). As expected, given the reduced recombination rate in subtelomeric regions, coldspots were particularly abundant in these regions. Specifically, subtelomeres contained 29.6% of the coldspots, although representing only 4.8% of the genome. This translates to a significant enrichment of coldspots in the 10% most terminal chromosomal regions compared to the proximal regions (Wilcoxon's W = 3; P = 5.62 × 10⁻⁴; Supplemental Fig. 2B). Permutation analyses revealed a significantly lower transposon density in coldspots, consistent across all classes (P-values < 0.01; Fig. 4A–D). Coldspots also showed a significant enrichment of genes related to transtransferase activity (P-values < 0.01, Supplemental Table 3), but no enrichment of specific sequence motifs (P > 0.05).

Figure 4.

Occurrence of transposable elements (TEs) in coldspots. The histograms show the density distributions of (A) DNA, (B) LTR, (C) SINE, (D) LINE, permuted means for 1283 random 30-kb genomic windows (number of, and average length of inferred coldspots). The red lines indicate the observed TE densities in the defined coldspot windows.

Open in new tabDownload PowerPointLink to figure

Association between recombination rate and base composition

We assessed the association between recombination rate and nucleotide composition at a 1-Mb scale to investigate the potential effects of gBGC in the wood white genome. The analysis showed that the recombination rate was overall not significantly correlated with the GC content (Spearman's rho = −0.07, P = 0.10; Fig. 5A). In contrast, this association was negative and marginally significant for the chromosomes involved in fusion/fission polymorphisms (Spearman's rho = −0.19, P-value = 0.04; Supplemental Table 4).

Figure 5.

Associations between the recombination rate and base composition (A) and gene density (B) in 1-Mb genomic windows. (A) Association between the recombination rate and the nucleotide composition (GC content). (B) Association between the recombination rate and the density of coding sequences (CDS).

Open in new tabDownload PowerPointLink to figure

Because we observed a significant reduction in recombination rate in subtelomeric regions, we investigated if those regions had deviating base composition. The analysis showed that there was a significantly higher GC content in subtelomeric regions (33.46%) compared to proximal regions (32.61%) (Wilcoxon's W = 1,109,810; P = 5.73 × 10⁻²²). The GC content was also slightly lower within hotspots (32.43%) and their 5 kb flanking regions (32.45%), compared to the genome-wide estimates (32.65%).

Associations between recombination rate and genomic features

We used multiple regression to assess the relative effect of six explanatory variables (GC content, gene density, DNA transposon, SINE, LINE, and LTR retrotransposon densities) on the recombination rate. This model revealed an overall significant association between the recombination rate and the explanatory variables (F(6) = 10.77, df = 603, P = 2.13 × 10⁻¹¹), but explained only a small fraction of the total variation (R² = 0.10, Adj R² = 0.09). Only SINE and LINE densities were significant explanatory variables (Table 1).

Table 1.

Summary of the results from a linear regression model with six explanatory variables included

Open in new tabDownload PowerPointLink to table

The recombination rate was neither associated with gene density in general (Spearman's rho = −0.04, P = 0.17; Fig. 5B), nor in chromosomes involved in segregating polymorphisms (P-value = 0.46, Supplemental Table 4). However, when analyzing different gene elements separately, we observed a significant reduction in recombination rate in exons (5.5 cM/Mb; Wilcoxon's W = 1.60 × 10¹⁰, P-value < 2.2 × 10⁻¹⁶) and 5′ UTRs (6.3 cM/Mb; Wilcoxon's W = 1.5869 × 10⁸, P-value < 2.2 × 10⁻¹⁶), and a significantly higher rate in introns (Wilcoxon's W = 7.3462 × 10¹², P-value < 2.2 × 10⁻¹⁶) compared to intergenic regions (7.5 cM/Mb) (Fig. 6).

Figure 6.

Recombination rate estimates in different gene elements and classes of transposable elements. The genome-wide recombination rate (7.37 cM/Mb) is indicated with the horizontal black dashed line. Orange bars indicate 95% confidence intervals.

Open in new tabDownload PowerPointLink to figure

The recombination rate was differently associated with different TE classes, both in strength and direction. DNA transposons (Spearman's rho = 0.09, P = 0.03, Fig. 7A) and SINEs (Spearman's rho = 0.29, P = 3.31 × 10⁻¹³, Fig. 7C) were positively associated, and LTRs (Spearman's rho = −0.11, P = 8.50 × 10⁻³, Fig. 7B) and LINEs (Spearman's rho = −0.19, P = 3.01 × 10⁻⁶, Fig. 7D) negatively associated with the recombination rate. In the chromosomes involved in segregating polymorphisms, only SINEs showed a significant correlation with the recombination rate (P-value = 3.96 × 10⁻⁵; Supplemental Table 4), and they were also the only significant explanatory variable in the linear regression analysis of these specific chromosomes (Supplemental Table 5).

Figure 7.

Associations between the recombination rate (1-Mb scale) and the proportion of four different TE classes; (A) DNA transposons, (B) LTRs, (C) SINEs, and (D) LINEs. The abundance of TEs was calculated as the fraction of each 1-Mb window occupied by each specific element.

Open in new tabDownload PowerPointLink to figure

Differences between TE classes were not restricted to the overall associations with recombination rate. The average recombination rate within each class of TEs also varied. The recombination rate was significantly lower for LINEs (6.6 cM/Mb; Wilcoxon's W = 5.233 × 10¹², P-value < 2.2 × 10⁻¹⁶) and DNA transposons (6.8 cM/Mb; Wilcoxon's W = 2.6104 × 10¹², P-value < 2.2 × 10⁻¹⁶), and significantly higher for SINEs (7.4 cM/Mb; Wilcoxon's W = 1.2212 × 10¹², P-value < 2.2 × 10⁻¹⁶) and LTRs (7.8 cM/Mb; Wilcoxon's W = 7.0468 × 10¹¹, P-value < 2.2 × 10⁻¹⁶), compared to the genome-wide estimates (Fig. 6).

Recent demographic trends are the consequence of post-Pleistocene population dynamics

The demographic trajectories were consistent across chromosomes, and historical fluctuations in N_e were in agreement with previous analyses, showing a gradual increase until approximately 10,000 generations ago (Talla et al. 2019). By using a method (SMC++) that accommodates unphased data from many individuals, we could infer trajectories in the recent past (Terhorst et al. 2017). We found a considerable decrease in N_e in the last 10,000 generations, a signal that Talla et al. (2019) could not pick up with PSMC (Li and Durbin 2011), which uses genomic information from a single individual (Mather et al. 2020). The recent decline in N_e is consistent with recurrent founder effects during the recolonization of Scandinavia after the LGP, and is compatible with the hypothesis of secondary contact and recurrent backcrossing underlying the observed chromosome number cline across Eurasia (Lukhtanov et al. 2018; Talla et al. 2019). The negligible differences in N_e obtained when using different sets of “distinguished” individuals reinforce the reliability of using SMC++ (Terhorst et al. 2017) for inferences on the modern demographic history of populations.

An important consideration regarding the time estimates of the demographic changes is that the assumption of univoltinity may be an oversimplification of the historically true scenario. Currently, wood whites are univoltine in the northern parts of the distribution range (Friberg et al. 2008), but voltinity gradually increases further south (Dantart 2003; Jeffcoate and Joy 2011). Historical changes in voltinity (e.g., Powell 1989; Margus and Lindström 2020) could affect the scale of inferences on population size evolution. However, the timing of the ice sheet retreat in Scandinavia is estimated at about 11,000 yr BP (Günther et al. 2018), matching the beginning of the observed decrease in N_e, indicating it is most likely a consequence of founder effects during the recolonization of Scandinavia from glacial refugia.

Genome-wide distribution of the recombination rate in wood whites

The genome-wide recombination rate varies considerably between different insects, from comparatively low in Diptera (<1 cM/Mb) (Beye et al. 2006) to exceptionally high in honeybees (19 cM/Mb), (Beye et al. 2006). We found that the genome-wide recombination rate in the Swedish wood white population was 7.37 cM/Mb. This is slightly higher than estimates from other lepidopterans like Bombyx mori (4.6 cM/Mb), Heliconius melpomene (5.5 cM/Mb), and Heliconius erato (6 cM/Mb) (Yasukochi 1998; Jiggins et al. 2005; Tobler et al. 2005), and substantially higher than in vertebrates, which range between 0.16 cM/Mb in the Atlantic trout to 3.17 cM/Mb in chicken (Beye et al. 2006). The recombination rate here is also somewhat higher than linkage map-based estimates for the Swedish population of L. sinapis (Näsvall et al. 2023). This is very likely a consequence of scaling with N_e inferred from demographic trajectories, rather than a true difference between the historical recombination rate (i.e., ρ) and direct estimates of r in families. Because we focus on recombination rate variation along chromosomes here, exact estimates of r are not crucial, but caution is of course needed when comparing our results to estimates of genome-wide recombination rate obtained using different approaches. The ZW sex-chromosome system in Lepidoptera, where males are homogametic, coupled with achiasmy (i.e., no recombination) in females (Suomalainen et al. 1973), could theoretically lead to a higher recombination rate on the Z Chromosome. The reason for that is that the Z Chromosome spends 2/3 of the time (on evolutionary timescales) in a male germline and 1/3 in a female germline (given equal sex-ratio). Our results do not support such an effect because no significant difference in recombination rate was observed between autosomes and the Z Chromosomes.

A clearly deviating recombination rate was found for Chromosome 16 (15.3 cM/Mb), and that explains why the association between chromosome size and recombination was nonsignificant. Negative associations between both variables have been observed in different taxonomic groups, for example, yeast (Kaback et al. 1992), humans and rodents (Jensen-Seaman et al. 2004), birds (Backström et al. 2010), cattle (Mouresan et al. 2019), and butterflies (Martin et al. 2019; Shipilina et al. 2022). Such a relationship is expected as crossovers are necessary for correct chromosome segregation during meiosis (Pardo-Manuel de Villena and Sapienza 2001; Smith and Nambiar 2020). Another genomic region of interest is the 1-Mb window on Chromosome 17 with the highest recombination rate in the entire genome (Fig. 2). Previous studies have shown that assembly errors may lead to biased estimates of LD and population recombination rate (Wall et al. 2022). However, both the female and male genome assembly in L. sinapis are collinear in that region (Höök et al. 2023), indicating that the high recombination observed in the vicinity is genuine.

Our analysis also showed that longer chromosomes tended to have a bimodal recombination landscape, with a reduced rate at the chromosome center and toward the chromosome ends. This pattern is consistent with findings in both monocentric organisms where crossing-over events are rare within the centromeres (Dapper and Payseur 2017), and holocentric species (e.g., Caenorhabditis elegans), where the recombination rate increases with the relative distance from the center of the chromosomes (Prachumwat et al. 2004). In holocentric lineages, the bimodal recombination landscape in larger chromosomes can be explained by the accommodation of multiple chiasmata. For example, observations in Psylla foersteri suggest that longer chromosomes can accommodate the formation of two simultaneous chiasmata, whereas shorter chromosomes cannot (Nokkala et al. 2004). In the cases where two chiasmata are formed in a single chromosome, crossover interference may prevent those from forming near each other and tend to drive them toward opposite ends of the chromosome (Otto and Payseur 2019). An additional, but not mutually exclusive explanation, is the formation of the “meiotic bouquet”—a stage in early meiosis characterized by the aggregation of chromosome ends close to the nuclear membrane, which can drive crossovers toward distal positions and reduce the recombination rate in the center of the larger chromosomes (Scherthan et al. 1996; Haenel et al. 2018).

Apart from a reduced recombination rate in the center of larger chromosomes, we also found a reduced recombination rate in terminal regions. This pattern was observed in almost all chromosomes (see exceptions below), irrespective of chromosome size and type. Note that the telomeres, which in Lepidoptera are 6–8 kb long tandem repeats of the motif (TTAGG)_n (Okazaki et al. 1993; Sasaki and Fujiwara 2000), were not assembled in the reference genome. Such decrease in the recombination rate in subtelomeric regions has previously been observed in Heliconius butterflies (Martin et al. 2019), as well as in yeast (DuBois et al. 2002; Barton et al. 2008) and flycatchers (Kawakami et al. 2017). A potential explanation for this pattern is that crossover initiation is prevented near the telomeres to minimize the risk for ectopic recombination between nonhomologous repeat sequences during meiosis (Smith and Nambiar 2020). As mentioned above, our results showed that the reduced recombination rate toward chromosome ends was not ubiquitous; eight chromosome ends (one for each of the following Chromosomes: 5, 6, 10, 11, 16, 25, 27, and 29) did not show a significant reduction in the recombination rate as compared to each respective intrachromosomal level. Four of these exceptions (one chromosome end for each of Chromosomes 5, 11, 25, and 27) coincide with recently identified fusion and fission polymorphisms segregating in the wood white Scandinavian population (Höök et al. 2023). These results show that fusion/fission events can have immediate effects on the distribution of crossovers within and between chromosomes.

Characterization of recombination hotspots and coldspots

The total number of hotspots (n = 3124) identified in the wood whites is similar to what was observed in Ficedula flycatchers using a comparable approach (Kawakami et al. 2017). The density of hotspots was, however, lower than in humans (n = 25,000–50,000 hotspots in an approximately five times larger genome) (Myers et al. 2005). Although the criteria for defining hotspots vary between studies, they all rely on contrasting the background and local recombination rates, making the results reasonably comparable. The distribution of hotspots is similar between wood whites and humans, as hotspots occurred mostly outside of genes (McVean et al. 2004). This is in contrast to birds, where genic and regulatory regions are enriched in hotspots (Singhal et al. 2015; Smeds et al. 2016; Kawakami et al. 2017). The frequency of recombination hotspots and coldspots was comparable in proximal regions, but the number of hotspots decreased, and coldspots were more frequent toward chromosome ends. A higher occurrence of recombination coldspots in terminal regions has previously been observed in yeast, which lack crossovers in chromosome ends altogether (Su et al. 2000; Barton et al. 2003). This observation is also in line with the observed decrease in average recombination rate in terminal regions in the wood whites and suggests that the recombination machinery is partly blocked from accessing the very ends of chromosomes.

We found that the recombination landscape in the wood white was highly variable. This is consistent with observations in humans (Mackiewicz et al. 2013), birds (Singhal et al. 2015), and dogs (Axelsson et al. 2012). However, other insects such as, for example, D. melanogaster (Comeron et al. 2012), for which detailed recombination maps are available, generally show less pronounced hotspots. Whereas the hotspot locations in humans are largely determined by the presence of sequence motifs associated with PRDM9 binding (Grey et al. 2011), little is known about what drives crossovers to occur at specific locations in organisms that lack a functional copy of PRDM9. To get preliminary information about potential mechanistic underpinnings of recombination rate variation in wood whites, we assessed if specific sequence motifs or gene categories were enriched in recombination hotspots and coldspots. The analyses did not reveal any associations for sequence motifs, and we found no enrichment of specific gene categories in hotspots. Coldspots, which experience almost no recombination, are interesting regions that can potentially be involved in the establishment of holocentromeres during meiosis in Lepidoptera. In these regions, we found a significant enrichment of genes with functions associated to transferase activity. Farnesyltranstransferase activity is a key step for the progression of cell division (Hussein and Taylor 2002) and correct attachment between the spindle and the kinetochores in humans (Moudgil et al. 2015). Although it is tempting to speculate that the active expression of farnesyltranstransferase might block the recombination machinery, their activity is expected to be transacting and they were found in only 13 of the coldspot regions. It is therefore not clear what the enrichment of farnesyltranstransferases in coldspots means from a mechanistic point of view, because other factors must underlie the repression of recombination in coldspots in general.

Associations between recombination, nucleotide composition, and gene content

The high-density recombination maps allowed us to investigate potential associations between the local recombination rate and different genomic features. Such information can be used to deduce the effects of recombination on base composition and/or regulatory mechanisms modulating the recombination landscape. A potential driver of an association between recombination and nucleotide composition is GC-biased gene conversion (gBGC), that is, the fixation bias favoring “strong” (G and C) over “weak” alleles (A and T) during recombination (Duret and Galtier 2009). This process mimics directional selection and can alter nucleotide composition between regions experiencing different recombination rates. We found no association between the overall recombination rate and nucleotide composition (GC content), and for chromosomes involved in fusion/fission polymorphisms there was a negative association. We also found that recombination hotspots had marginally lower GC content (the opposite is expected with gBGC; Kawakami et al. 2017), consistent with a limited effect of gBGC in butterflies (Boman et al. 2021), and in contrast to findings in several other systems like humans (Fullerton et al. 2001; Meunier and Duret 2004), mice (Clément and Arndt 2013), flycatchers (Kawakami et al. 2017) and fruit flies (Marais et al. 2001), as well as plants (Muyle et al. 2011) and yeast (Gerton et al. 2000; Kiktev et al. 2018). In yeast, plants, and humans, mononucleotide runs have been proposed as drivers of the association between recombination rate and nucleotide composition, counteracting the effects of gBGC on long-term nucleotide composition (Bagshaw et al. 2008; Leushkin and Bazykin 2013; Heissl et al. 2019). However, the enrichment analysis showed no overrepresentation of mononucleotide runs in either hotspots or coldspots, indicating that other factors mediate the recombination landscape in butterflies.

In many organisms, for example, mice (Paigen et al. 2008) and different plant species (Gaut et al. 2007; Tiley and Burleigh 2015), recombination occurs more frequently in gene-dense regions, but we did not find such an association in the wood whites. However, our data showed that the recombination rate was significantly reduced in exons and 5′ UTRs compared to introns and intergenic regions. This is in line with findings in other insects (Wallberg et al. 2015; Jones et al. 2019), and humans (McVean et al. 2004), where recombination hotspots concentrate in the vicinity of, but not within, coding and regulatory regions. The slightly elevated recombination rate in introns is consistent with findings in the holocentric nematode C. elegans (Prachumwat et al. 2004), but opposed to observations in Drosophila (Carvalho and Clark 1999) and humans (Comeron and Kreitman 2000). Taken together, this indicates that recombination occurs within genes in butterflies, but that crossovers are partly inhibited in coding sequences which might lead to a slightly elevated rate in introns.

Different TE classes show contrasting association patterns with the recombination rate

Associations between recombination and TE densities have mainly been investigated in organisms with defined centromeres (Kent et al. 2017), although data from holocentric species are scarce (but see Lavoie et al. 2013; Baril and Hayward 2022; Smolander et al. 2022). To assess the relationship between the recombination rate and genomic features in the wood white, we used TE annotation information previously developed for the species (Höök et al. 2023). This analysis revealed differential associations between the recombination rate and the abundance of distinct TE classes. DNA transposons and SINEs were positively associated, although LINEs and LTRs were negatively associated with the recombination rate. Under the assumption that TE insertions in general are slightly deleterious, we would expect a negative correlation between the recombination rate and TE abundance (Kent et al. 2017) as a consequence of more efficient purging of deleterious insertions in regions of high recombination rate (Bartolomé et al. 2002; Wright et al. 2003). However, given that recombination is initiated by a double-strand break, it is possible that certain TEs are used as a template for the repairing process, increasing their frequency in regions of higher recombination rates (Onozawa et al. 2014). SINEs have, for example, been shown to use DNA breaks to integrate back into the genome after replication (Singer 1982). Associations between TEs and the recombination rate are hence expected to depend on the TE class and the focal study system. For example, Alu elements (a subfamily of SINEs) in humans accumulate in regions with elevated recombination (Witherspoon et al. 2009), and SINEs are positively associated with the recombination rate in the painted lady (Vanessa cardui) (Shipilina et al. 2022). Similarly, DNA transposons are enriched in regions of high recombination rate in C. elegans (Duret et al. 2000). The causality of such associations is not easy to establish. In cases in which TE proliferation has deleterious fitness effects, we expect a negative association between TE abundance and recombination rate. However, Alu elements in humans have been shown to be drivers of an increase in the local recombination rate, possibly caused by Alu elements mimicking the action of short recombinogenic motifs (Witherspoon et al. 2009). It is also possible that other underlying factors affect the TE and recombination rate distributions similarly; both the TE proliferation and the recombination initiation machinery, for example, target open chromatin more easily (Kawakami et al. 2017). In the wood whites, the average recombination rate estimates within TE classes did not deviate considerably from the genome-wide average. However, these comparatively minor differences can indicate differences in the selective pressure against insertion of specific families of TEs. This does not seem to be related with the length of the TEs, as SINEs considerably shorter than LINEs and LTR elements showed an average recombination rate between the longer types. For chromosomes involved in fusion/fission polymorphisms, the only association between recombination rate and TE density was found for SINEs (positive). The lack of a negative association with LINEs and LTRs in the rearranged chromosomes may potentially be a consequence of reduced efficiency of selection (Dolgin and Charlesworth 2008), but could also just reflect reduced power in the analysis with the smaller set of data points.

In summary, the different TE classes showed different associations with the recombination rate. This is consistent with findings in other studies and may point toward similar determinants in holocentric organisms compared to those with defined centromeres. However, the causality needs further study through detailed characterization of crossover regions in large pedigrees.

Genome assembly

The wood white genome assembly used as reference was developed by Höök et al. (2023). In brief, one mated adult female wood white was caught in Sweden and reared for egg laying. One male offspring pupa was sampled and flash frozen in liquid nitrogen. The sample was divided to create a 10x Genomics Chromium Genome-library and a Dovetail Hi-C library. For 10x sequencing, DNA was extracted using a modified HMW salt extraction method (Aljanabi and Martinez 1997). Tissue for Hi-C sequencing was disrupted in liquid nitrogen. Library preparation, sequencing, and genome assembly were performed by NGI Stockholm. Sequencing was performed on Illumina NovaSeq 6000 with 2 × 151 setup. 10x linked-reads were assembled with 10x Genomics Supernova v2.1.0 (Weisenfeld et al. 2017). Hi-C reads were processed with Juicer v1.6 (Durand et al. 2016a) and used for scaffolding the 10x assembly with 3DDNA v.180922 (Dudchenko et al. 2017). Resulting assemblies were reviewed with Juicebox v1.11.08 (Durand et al. 2016b). In addition to minor corrections to the initial assembly, two chromosome-sized scaffolds were merged. The assembly was finalized with the script “run-asm-pipeline-post-review.sh” from 3DDNA v.180922 (Dudchenko et al. 2017).

Gene annotation

Gene annotation liftOver was performed by aligning wood white protein queries generated by Talla et al. (2017) to the current version of reference genome, using spaln 2.4.0 (Iwata and Gotoh 2012) with parameters -Q7 -LS -O7 -S3.

TE annotation

Repetitive consensus sequences were predicted de novo using RepeatModeler 1.0.11 (Bao et al. 2015). TEs characterized as unknown were submitted to CENSOR (Bao et al. 2015) for annotation, where hits with a score < 200 were removed. Predicted sequences were matched against gene annotations using DIAMOND BLAST 2.0.4 (Buchfink et al. 2021), to correct for annotation errors (bitscore > 100). TEs were annotated in the genome with RepeatMasker 4.1.0, using the predicted library of consensus sequences in wood white and previously characterized TEs in Heliconius melpomene in the RepeatMasker library 4.0.8 (Bao et al. 2015).

Sampling of individuals

Adult male wood whites were collected across Sweden during June and July 2020 (Supplemental Table 6). Sex was determined in situ based on two sexually dimorphic characters: the presence of a black apical forewing spot and white coloration of the ventral part of the antennae in males. Sampled individuals were preserved in ethanol and frozen at −20°C.

DNA extraction

DNA was extracted following two protocols. In both, the dissected tissue was digested overnight in Laird's buffer and homogenized with 20 μL of Proteinase K (20 mg/mL, >600 mAU/mL), and incubated with RNase A at 37°C for 30 min. DNA from thoraces was extracted using salt extraction; 300 μL of 5M NaCl was added, followed by centrifugation for 15 min at 13,000 revolutions per minute (rpm). Three washing steps were completed with one volume of 70% ethanol and centrifuging for 5 min at maximum speed. The pellet was air-dried and then resuspended in 30 μL of Milli-Q H₂O. For the abdomens, a phenol-chloroform extraction protocol was used. Two cycles of phenol:chloroform:isoamyl alcohol (25:24:1) addition and centrifugation for 5 min at 13,000 rpm were completed, plus a third cleaning cycle using only chloroform. DNA precipitation was achieved by adding 2× volumes isopropanol + 0.1× 3M NaAc, incubating at −18°C overnight and centrifuging 15 min at 13,000 rpm. The final pellet was resuspended with 30 μL of Milli-Q H₂O. DNA purity was assessed with NanoDrop, and concentration measured with Qubit DNA Broad Range.

Sequencing

84 individuals with the highest DNA quality representing the distribution range were selected. Library preparation for all 84 samples using the TruSeq PCR-free kit followed by multiplexing, and sequencing on two NovaSeq 6000 S4 lanes with 2 × 150 bp reads, were performed by NGI Stockholm.

Read trimming

Illumina adapters were trimmed by eliminating 15 bp on each end of the raw reads with cutadapt 1.9.1 (Martin 2011), filtered on Q-score < 30 and minimum length of 30 bp. Read quality postcleaning was assessed with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Before filtering, 4.3 million reads per sample were obtained on average, and 2.5% were filtered out.

Mapping and filtering

Paired-end reads were mapped to the reference genome with BWA v0.7.17 (Li and Durbin 2009). SAMtools v1.10 (Li et al. 2009) was used to select reads with paired information. Duplicates were eliminated with MarkDuplicatesSpark as implemented in GATK v4.1.4.1 (McKenna et al. 2010) with the ‐‐remove-sequencing-duplicates option.

Variant calling and filtering

GATK's v4.1.4.1 HaplotypeCaller (McKenna et al. 2010) was used for variant calling. All chromosomes and individuals were called independently, and the resulting 84 files for each chromosome were grouped into a single VCF file with the GATK's v4.1.4.1 Combine_gVCF and Genotype_gVCF (McKenna et al. 2010). Total variant count was obtained with the stats option in BCFtools v1.10 (Li 2011). Variants were filtered to have a minimum minor allele count of two, per site depth between 10 and 50, minimum per site quality of 30, and <20% per site missing data, using VCFtools 0.1.15 (Danecek et al. 2011). Insertions and deletions were removed with the ‐‐remove-indels option. The number of remaining sites in each chromosome was counted again with the stats option in BCFtools v1.10 (Li 2011). Initial variant calling resulted in 51,189,479 markers (11,055,543 indels + 40,133,936 SNPs), 10,565,404 SNPs remaining after filtering. This represents a genome-wide average of ∼17.6 SNPs/kb, given the ∼ 600-per Mb total length of the reference genome.

Inference of the demographic history

The demographic history for each chromosome was inferred with SMC++ (Terhorst et al. 2017) using six “distinguished” individuals. These were selected based on high sequencing coverage after filtering and to be representative of the geographic distribution. SMC++ vcf2smc was run for each “distinguished” individual, taking the allelic state of all other samples into account for the final output. SMC++’s estimate was used to infer the demographic history of each chromosome by incorporating the genealogical information of all six runs and computing the composite likelihood. The per base, per generation mutation rate was set at 2.9 × 10⁻⁹, based on mutation frequency estimations in H. melpomene (Keightley et al. 2015). Invariable regions such as centromeres must be masked before demographic inferences, as they can interfere with the signal. Because Lepidoptera are holocentric and lack defined centromeres, a cutoff value of 150 kb was set instead, so that any longer invariable region was considered as missing data and discarded for the demographic inference. The demographic trajectories were inferred for the last 5 × 10⁶ generations. To test the accuracy of the inferences, the demographic analyses were performed again using a randomly selected set of six distinguished individuals. The most recent estimates of N_e, which are key for scaling the population recombination rate with pyrho (Spence and Song 2019), were used to assess if the choice of distinguished individuals significantly changed the output.

Recombination rate estimation

The chromosome-specific demographic trajectories were used together with the VCF files to obtain high-resolution recombination maps using pyrho (Spence and Song 2019). An algorithm implemented in LDpop (Kamm et al. 2016) was used to compute two-loci likelihoods under the coalescent with recombination using the chromosome-specific demographic trajectories as input. The same mutation rate as before was used, together with the parameter values n = 168 (twice the number of diploid individuals), and N = 210 (25% larger than n, as recommended in the manual). A relative tolerance (‐‐decimate_rel_tol) of 0.1 was used with the ‐‐approx flag, recommended for large data sets. Different window sizes (maximum distance between SNPs) and block penalties (determinant of the smoothness of the curve) were tested for each chromosome, and those with the highest correlation between the likelihood tables and simulated data at different scales (1 bp, 10 kb, and 100 kb) were selected. For all chromosomes, the best block_penalty was 25, and the best suitable window_size ranged between 50 and 100. The look-up table and the final VCF file for each chromosome were used to infer the local recombination rate using the most appropriate parameter values and the ‐‐fast_missing flag. The per base, per generation recombination rate between each pair of SNP markers was obtained in the end. Our short-read data did not allow for the elimination of heterozygous individuals for the identified rearrangements. However, fusion/fission polymorphisms in the species are comparatively recent with respect to the buildup of LD patterns, hence they should not be a major cause of biases in the inferences and scaling from ρ to r. A positive association between the recombination rate and marker density was observed at a 1-Mb scale (Supplemental Information; Supplemental Fig. 3).

Distribution of recombination rate variation and identification of hotspots and coldspots

Regional recombination rate estimates were obtained in two window sizes (100 kb and 1 Mb) by calculating the weighted average recombination rate between markers, accounting for the distance between them. We used a simulation-based approach to establish thresholds for hotspot identification. Coalescent simulations were used to assess how much of the variation inferred by pyrho was caused by our imputation strategy. msprime 1.1.1 (Kelleher et al. 2016) was used to simulate 99 independent sequences, each 100 kb long, and the VCF files for each sequence were concatenated. The parameters for the simulations included a flat recombination landscape (7.37 cM/Mb, the genome-wide average) and the same mutation rate used for SMC++, together with an approximation of the demographic trajectory inferred from the empirical data; exponential growth in the interval 10⁴–10⁶ generations BP, exponential decline 10²–10⁴ generations BP, and stable size in the last 10² generations. At each time point, N_e was calculated as the average of the estimates across chromosomes. The simulated genomic sequences were analyzed in pyrho according to the steps described earlier. The resulting regional recombination rates predominantly oscillated in the range 2–8 cM/Mb, with a maximum of 20 cM/Mb, roughly threefold higher than the genome-wide average (Supplemental Fig. 4). Recombination peaks occurring at boundaries between independently obtained VCF files were omitted. Therefore, simulations allowed us to establish a lower threshold for hotspot consideration. For the final (more conservative) criteria, we defined hotspots as regions with a recombination rate higher than 25 cM/Mb, between 750 and 10,000 bp long, and showing a 10-fold increase over the regional “background” recombination rate. This “background” recombination rate was established for each 100-kb window, as the average rate for the focal window and the two flanking windows, in total 300 kb. To avoid biases resulting from erroneously called variants, hotspots that included less than four markers (i.e., three intervals) were discarded. Recombination coldspots were identified as regions with a recombination rate 10-fold lower than the genome-wide average, and including ≥4 markers with no length limitations.

To assess the distribution of recombination hotspots and coldspots, we defined the subtelomeric regions as the last 500 kb on each chromosome end. The rest of the genome was considered proximal. Flanks of recombination hotspots were defined as the 5-kb segments at each side of the defined hotspots.

Associations between recombination rate with GC content, gene density, and classes

A linear regression model was applied to disentangle the explanatory potential of each genomic feature for the recombination rate variation at a 1-Mb scale, using the lm function in base R (R Core Team 2021). The model had recombination rate as the response variable, and the explanatory variables were GC content, gene density, and the relative abundance of four TE classes (DNA transposons, and LINE, SINE, and LTR retrotransposons). Potential associations between local recombination rate and the different genomic features were also analyzed with Spearman's rank correlation tests using the corr option in base R (R Core Team 2021).

BEDTools 2.29.2 (Quinlan and Hall 2010) maskfasta was used to select all annotated exons and TEs. Base composition for specific regions (corrected for masked positions) was obtained with BEDTools nuc option. Gene and TE densities were calculated as the proportion of a region covered by the annotated sequences of each category.

To assess the variation in the recombination rate within genomic features, we estimated the recombination rate within each TE class, different regions in protein-coding genes (exons, introns, and upstream regions) and intergenic sequences. Because 5′ UTRs were not annotated, we used the 100-bp upstream of each gene's first exon, a conservative selection to represent 5′ UTRs (Chen et al. 2011, 2014). To avoid biases caused by differential selective pressures, only the intervals between SNPs within a feature were considered; for example, for a given exon, if an interval between markers overlapped completely with the predicted exon, it was retained; if the overlap was only partial (including part of the exon sequence but also part of a neighboring intron or UTR), it was discarded. This overlap was checked with BEDTools 2.29.2 (Quinlan and Hall 2010) intersect option using the flag -f 1.

Permutation test to assess TE density in the coldspots and hotspots

We ran 10,000 permutation tests in R v4.2.1 (R Core Team 2021) with windows reflecting the count and average length of recombination hotspots and coldspots to obtain neutral distributions of TE abundances. We then compared these distributions with the empirical data and assessed the significance of our findings with two-sided tests.

Gene Set Enrichment Analysis (GSEA) for the genes present in the coldspot and hotspot regions

We assessed potential enrichment of functional categories in the hotspot and coldspot regions using topGO v2.44.0 (https://bioconductor.org/packages/release/bioc/html/topGO.html) in R v4.2.1 (R Core Team 2021). We used the annotated gene set with Gene Ontology (GO) terms associated to molecular function and Fisher's exact test to test for significance, with the default algorithm (“weight01”) to account for the hierarchical structure of the GO terms (Alexa et al. 2006). We adjusted the P-values with Benjamini–Hochberg's method of multiple test correction (p.adjust(x, method = “fdr”)). HOMER v4.11 (Heinz et al. 2010) was used to check for enriched sequence motifs within the coldspots and hotspots by comparing the abundance of oligonucleotide sequences with background genome-wide frequencies, with autonormalization of up to 3-mers.