Local adaptation of life cycles in a butterfly is associated with variation in several circadian clock genes

Many insects exhibit geographical variation in voltinism, the number of generations produced per year. This includes high‐latitude species in previously glaciated areas, meaning that divergent selection on life cycle traits has taken place during or shortly after recent colonization. Here, we use a population genomics approach to compare a set of nine Scandinavian populations of the butterfly Pararge aegeria that differ in life cycle traits (diapause thresholds and voltinism) along both north–south and east–west clines. Using a de novo‐assembled genome, we reconstruct colonization histories and demographic relationships. Based on the inferred population structure, we then scan the genome for candidate loci showing signs of divergent selection potentially associated with population differences in life cycle traits. The identified candidate genes include a number of components of the insect circadian clock (timeless, timeless2, period, cryptochrome and clockwork orange). Most notably, the gene timeless, which has previously been experimentally linked to life cycle regulation in P. aegeria, is here found to contain a novel 97‐amino acid deletion unique to, and fixed in, a single population. These results add to a growing body of research framing circadian gene variation as a potential mechanism for generating local adaptation of life cycles.


| INTRODUC TI ON
Organisms in temperate regions often experience a radically different environment depending on the time of year, and face a major evolutionary challenge in adapting their life cycles to this seasonality. Strong seasonality will favour genotypes that execute particular life cycle events (e.g., mating; migration; entering or exiting dormancy) at the appropriate times of year. Because the seasonal cycle is not the same everywhere, but varies with local climate, optimal life cycles also vary geographically. Hence, seasonality provides an opportunity to study divergent adaptation across populations (Bradshaw et al., 2004;Conover & Schultz, 1995;Posledovich et al., 2015;Stinchcombe et al., 2004).
In insects, a central life history trait is voltinism, the number of generations produced by a population in a year. Typically, in a population with one generation per year (univoltinism), every individual at some (species-specific) stage of life enters diapause, a state of hormonally suspended development and metabolic suppression, in anticipation of winter (Koštál, 2006;Tauber et al., 1986). In contrast, a sufficiently long warm season may permit one or several additional, nondiapausing generations to occur each year (bivoltinism or multivoltinism, respectively). Producing a single versus multiple generations per year is a nontrivial difference, as it is associated with drastic differences in generation time (Roff, 1980), and may necessitate or drive associated adaptations in terms of regulation of diapause induction (Bean et al., 2012;Lindestad et al., 2019), life history traits such as development rate and body size (Fischer & Fiedler, 2002;Masaki, 1972;Mousseau & Roff, 1989), and developmental pathway expression (Aalberg Haugen & Gotthard, 2015;Kivelä et al., 2013). Nonetheless, many insect species differ geographically in voltinism, typically in relation to the local length of the warm season (Braune et al., 2008;Faccoli & Bernardinelli, 2014;Hart et al., 1997;Hu et al., 2012;Levy et al., 2015).
Expression of a particular local voltinism is often tied to photoperiodism; that is, how the population responds to differences or changes in daylength. The most widely described regulatory mechanism underlying differences in voltinism is local adaptation of the critical daylength (CD) below which diapause is induced (Bean et al., 2012;Grevstad & Coop, 2015;Lindestad et al., 2019). CD tends to be locally adapted: in particular, a positive correlation with latitude has been recorded across a diverse set of taxa (Hut et al., 2013).
Generally speaking, a lower CD means that diapause will be induced later in the year, extending the insect's active season, which allows additional generations to take place in the same year. However, because voltinism expression depends on the interplay of several environmental and physiological factors (Grevstad & Coop, 2015), there is not a simple correspondence between a population's CD and its voltinism.
Genetic studies have shown that interpopulation differences in diapause induction are generated by a combination of variation at loci of small and large effects (Ragland et al., 2019). In several cases (Dalla Benetta et al., 2019;Levy et al., 2015;Pruisscher et al., 2018;Tauber et al., 2007;Yamada & Yamamoto, 2011), the detected large-effect loci are components of the circadian clock, which exhibit a widely recognized, but controversial and poorly understood, link with photoperiodism across many insect systems (Koštál, 2011;Meuti & Denlinger, 2013). To determine how widespread this association between life cycle adaptation and circadian gene variation is, there is a need for broad analyses that survey genome-wide variation potentially associated with these traits. Such analyses will also be aided by a better understanding of the historical and nonselective genetic backgrounds against which selection on life cycle traits has been acting.
The present study explores genetic variation in the speckled wood (Pararge aegeria), a woodland-associated satyrine butterfly that varies greatly in voltinism across its pan-European range. The northernmost P. aegeria populations are univoltine; southern populations have many overlapping generations per year and appear not to diapause at all (Nylin et al., 1995). Common-garden studies have demonstrated that this life cycle variation is to a large extent generated by local adaptation of photoperiodic plasticity (Lindestad et al., 2019;Nylin et al., 1995). In mainland Scandinavia, P. aegeria populations shift from univoltine in central Sweden to bivoltine in southern Sweden and Denmark (Figure 1a). The southern Swedish populations appear to be recently established, as no record of them exists before the 1930s, and they are still separated from the central Swedish populations by a distribution gap wherein P. aegeria is rare or absent (Eliasson et al., 2005;Nordström, 1955). Historical modelling from microsatellite data shows that both central Swedish and south Swedish populations probably immigrated from the south via Denmark (Tison et al., 2014), which would mean that south Swedish populations went extinct or nearly extinct after the initial northward colonization, and later rebounded either from local stock or through a second migration event via Denmark. In addition, bivoltine populations exist on the large Baltic islands of Öland and Gotland, on the same latitude (and hence experiencing the same photoperiod) as some of the univoltine mainland populations. The Öland population is genetically depauperate, while the Gotland population maintains quite high genetic diversity (Tison et al., 2014). The colonization histories of both island populations are unknown.
Previously, a study of genome-wide differentiation patterns between a univoltine population (Sundsvall) and a bivoltine population (North Skåne) at either end of the range of P. aegeria in Sweden (Pruisscher et al., 2018) identified 15 genomic regions showing signs of divergent selection. Single nucleotide polymorphism (SNP) genotypes in some of these regions-including nonsynonymous substitutions in two circadian loci, period and timeless-also correlated with diapause incidence in a laboratory test of interpopulation hybrids, suggesting that variation at these candidate loci contributes to maintaining voltinism variation across the studied geographical region by affecting photoperiodic plasticity. Furthermore, genotyping of three additional populations showed that the cline in voltinism from north to south is accompanied by a gradient of allele frequencies at the identified candidate loci.
However, these previous results cannot explain all of the life cycle variation seen in the studied populations. The clearest indication of this is the Gotland island population, which is bivoltine and exhibits a relatively low CD (Aalberg Haugen & Gotthard, 2015;Lindestad et al., 2019), despite closely resembling the univoltine, high-CD mainland populations in terms of variation at the 15 candidate loci (Pruisscher et al., 2018). In other words, if there is indeed a causal link between the identified candidate SNPs and diapause incidence in the studied populations, there must in addition be a significant amount of genetic variation modifying this trait that remains to be identified. Extending the genome-wide search to include more intermediate-latitude populations, especially populations with contrasting voltinism and different colonization histories, should help identify such variation.
Here we build on these earlier results in several ways. First, we assemble a high-quality genome for P. aegeria in order to improve our analyses. Second, we use whole-genome resequencing to reconstruct phylogeographical relationships between nine Scandinavian P. aegeria populations, and infer historical gene flow events. Third, we scan the populations for genomic regions showing signs of divergent selection potentially related to variation in life cycles, across both north-south and island-mainland voltinism clines. Our findings contribute to a more complete picture of the selective and historical processes underlying variation in insect life cycles, including a potential contribution from novel variation in circadian genes that may help resolve previous incongruities in this study system.

| Genome assembly and annotation
High-molecular-weight DNA was extracted from a single Pararge aegeria female from Skåne in southern Sweden, and used for constructing a 10× Genomics Chromium Genome library. Library preparation, sequencing and genome assembly were performed at the National Genomics Infrastructure, SciLifeLab, Stockholm, Sweden.
Genome assembly was conducted using supernova 1.2 with the phase option enabled; 12 separate assemblies were carried out, each using a different proportion of the 10× sequencing data (from 15% to 100%). Each generated assembly was assessed for two metrics: basic assembly stats, using quast 4.0 (Gurevich et al., 2013), and genic content, using busco version 3 (Waterhouse et al., 2018) with the database eukaryota_odb9. The genome assembly deemed optimal according to these metrics was then scaffolded using a 5-kb insert mate-pair library generated from the same P. aegeria individual, after quality filtering and mate pair-specific filtering using nextclip version 1.3 (Leggett et al., 2014), using the scaffolding software besst version 2.0 (Sahlin et al., 2014). The genome was polished using genomic data from a single individual via the software pilon version 1 (Walker et al., 2014). Next, the genome was collapsed to a haploid copy using haplomerger2 version 20180603 (Huang et al., 2017), and again assessed for genic content using busco, this time with the database insecta_odb9. Then, to improve gene model regions, the F I G U R E 1 Geographical background and overall patterns of genetic diversity across populations. Except where noted otherwise, northern cluster populations are written in red, southern cluster populations in blue, and the two island populations in orange and brown, respectively. (a) Map of Pararge aegeria distribution across Sweden and Denmark. Light-shaded and dark-shaded regions represent predominantly univoltine and bivoltine populations, respectively (adapted from Lindestad et al., 2019). Points mark sampling locations; point size reflects sample size. Population numbers show the critical daylength for diapause induction, in hours (see text for details). (b) Phylogeographical relationships as inferred by treemix. Branch lengths represent the amount of genetic drift. Dotted arrows represent inferred migration events; arrow thickness represents the proportion of the recipient population's ancestry that derives from the donor population. (c) Average genome-wide nucleotide diversity (π) and Tajima's D per population. (d) Average genome-wide F ST between all pairs of populations. Point shape and colour reflects comparison type: large blue circles, southern mainland populations; red squares, northern mainland populations; purple triangles, northern vs. southern populations; small orange circles, Öland versus all other populations; brown diamonds, Gotland vs. mainland populations assembly was super-scaffolded using the proteome from the Bicyclus anynana assembly by Nowell et al. (2017). This proteome (Bicyclus_ anynana_nBa.0.1._-_proteins.fa.gz) was downloaded from Lepbase (Challis et al., 2016). Superscaffolding was carried out using mespa (Neethiraj et al., 2017), which is a wrapper for the high-performing protein-to-genome aligner spaln2 (Iwata & Gotoh, 2012). Finally, the genome was annotated using the B. anynana proteome, as well as a pre-existing RNA sequencing (RNAseq) data set from P. aegeria, by running braker2 (Hoff et al., 2016;Stanke et al., 2006Stanke et al., , 2008 on an assembly that had been soft-masked for repetitive content using red (Girgis, 2015). Functional annotation was generated using the eggnogversion 5.0 online server with default settings (Huerta-Cepas et al., 2019).

| Pooled resequencing of populations
DNA for pool-sequencing was extracted from 201 P. aegeria individuals sampled from eight populations across Sweden and Denmark ( Figure 1a; Table 1). Most individuals (n = 168) were wild adults caught in 2014-2016; of these, 12 female butterflies could not be used, but as each female was already mated when caught, one offspring per female was used instead. In order to improve sample size for two of the populations, Öland and Stockholm, 33 individuals were added that had been collected in 2010-2011 for a previous study (Tison et al., 2014). Both sexes were represented for all populations. The tissue used was either half of the thorax (for adults) or the posterior end of the abdomen (for pupae).
Extraction was carried out separately per individual, using the standard "cell and tissue DNA" kit on a KingFisher Duo Prime purifier (ThermoFisher Scientific), with added RNAse A to avoid RNA contamination. DNA concentration was measured using a Qubit 2.0 flourometer (ThermoFisher Scientific), and samples were pooled per population using an equal amount of DNA from each individual, yielding 4-4.5 µg DNA per population pool. Pooled samples were run on a 2% agarose gel stained with GelRed to visually ascertain that DNA fragmentation was minimal. Library preparation and sequencing (Illumina HiSeq) was performed by SciLifeLab (Uppsala, Sweden), using 150-bp paired-end reads with 350-bp insert size.
In addition to the eight population read sets generated here, an additional P. aegeria read set (100-bp paired-end reads; 450-bp insert size) from Sundsvall in northern Sweden was used in our analyses.
The Sundsvall data set was generated from 22 adult P. aegeria individuals caught in 2011 for a previous study, using similar extraction, pooling and sequencing methods (Pruisscher et al., 2018).

| Data set preparation
Raw sequencing files were filtered for polymerase chain reaction (PCR) duplicates using the clone_filter script of stacks 1.21 (Catchen et al., 2013), and Illumina sequencing adapters were removed using bbduk2 (Bushnell, 2015). All nine PoolSeq read sets were mapped to the P. aegeria genome using nextgenmap version 0.4.10 (Sedlazeck et al., 2013) at 90% identity. Alignments were filtered using samtools version 1.6 (Li et al., 2009), keeping only properly paired reads with a map quality of 20 or higher. Read depth (RD) exceeded pool size over large parts of the alignments, leading to a risk of unequal read sampling of sequenced individuals, and hence skewed estimates of nucleotide diversity (Hoban et al., 2016). Because the analyses relied on accurate diversity estimates, and for consistency across analyses, we applied a subsampling approach: reads were sampled without replacement up to a population-specific target RD. This meant that genomic positions where raw RD was lower than the target RD were discarded, and so, raw RD distributions were examined in order to pick an appropriate TA B L E 1 DNA sample details. Populations considered as belonging to the northern or southern phylogeographical clusters are marked (N) or (S), respectively. Pool size is the total number of individuals used; of these, some (italicized numbers in parentheses) were caught for a previous study (Tison et al., 2014). For the Sundsvall pool, DNA extraction and sequencing had also been done previously (Pruisscher et al., 2018). For populations represented by two sampling locations, coordinates for the site contributing the largest part of the sample are listed first. Sex is the number of individuals of each sex used (this information was missing for a few individuals, marked "?"). DNA used refers to the total mass of DNA of the pooled sample target RD that allowed a majority of the data to be kept while still ensuring accurate allele frequency estimates. The chosen cutoff was the 5th percentile of raw RD for each respective population, or a minimum of 20.
Before subsampling, each alignment was converted into a population-wise pileup file using samtools. To avoid spurious SNP calling, indels were filtered out (with a 5-bp margin on each side) using scripts from the popoolation version 1.2.2 package (Kofler, Orozco-terWengel, et al., 2011). Next, also using popoolation, reads were subsampled from each pileup up to the target RD. Pileup positions in the 99th percentile of per-population RD were also excluded.
The resulting nine subsampled population-wise pileups (details in Appendix S1: Table S1) served as the input for all downstream analyses. A PHRED quality score cutoff of 20 and a minimum count of 2 for SNP alleles were used in all subsequent analyses.

| Population history and connectivity analyses
To estimate overall population differentiation, the average genomewide fixation index (F ST ) was calculated for all possible population pairs using popoolation2 version 1.201 . All nine input pileups were joined together as separate columns in a single mpileup file, which was then converted into the sync format. SNPs were called on this sync file, yielding about 7.8 million SNPs (the precise number of SNPs called in each analysis is shown in Appendix S1: Table S2). F ST was calculated on the called SNP variation in nonoverlapping windows across the whole genome. The window size used was 5000 bp, as this was found to be large enough to prevent noisy results, but small enough to still capture population differences.
The mean F ST across all windows for each pairwise population comparison was then calculated. Here and in all other analyses, F ST was calculated from allele frequencies as per Hartl and Clark (2007), as implemented by default in popoolation. Additionally, overall genetic diversity in each population was quantified by measuring nucleotide diversity (π) and Tajima's D across each single-population pileup, again in 5,000-bp windows, using popoolation. The D statistic measures the relationship between nucleotide diversity and the number of variable sites, and can provide information on demographic history and selection dynamics (Tajima, 1989a(Tajima, , 1989b. For both the F ST , π and Tajima's D analyses, only genomic windows where at least 50% of positions met the RD requirements were analysed. In order to investigate population history and potential gene flow, a phylogeographical analysis was conducted using treemix version 1.13 (Pickrell & Pritchard, 2012). treemix takes genome-wide allele frequency distributions as input, reconstructs a bifurcating phylogeny and identifies population pairs that share more allele frequency variation than the reconstructed tree can account for. The software then allows these incongruences to be resolved by sequentially adding inferred migration events between points on the tree, where each migration event is modelled as a unidirectional donation of a specific percentage of a population's allelic variation. Input allele frequencies were extracted from the same nine-population sync file as for the F ST analysis. treemix was run five times, each time allowing for one additional inferred migration event, from zero migration events to four. Because previous results indicate that all Swedish regions were settled from the south (Tison et al., 2014), Copenhagen was set as the root population in all runs.

| Scanning for footprints of selection (F ST and π)
As a first step towards identifying genomic regions potentially in- Lepidoptera (moths and butterflies) as the focal taxon.

| Scanning for loci associated with CD
To further identify genomic regions potentially associated with voltinism variation, while also accounting for the underlying demographic relationships among the studied populations, we applied a genotype-phenotype association approach using baypass (Gautier, 2015). For this analysis we focused specifically on CD for pupal dia- . On the assumption that CD is a linear function of temperature within the studied range, population means were calculated from a regression model with each estimate as one data point, population as a factor and temperature as a covariate (Appendix S1: Figure S1). The fitted per-population CDs at 18°C (Figure 1a) were then used as the phenotypic trait values for the baypass analyses.
baypass uses a Bayesian approach to identify loci whose allele frequencies correlate with a phenotypic trait across populations, while controlling for the underlying population covariance structure. The input used was the same allele frequency table as for the treemix analysis. (Although the treemix and baypass analyses were run independently, the population structure inferred by baypass and used as the basis for its analysis was highly similar to the one produced by treemix.) CD phenotype values were scaled using the -scalecov command, and baypass was run using the standard (default) model.
The analysis was run in triplicate, and for each analysed SNP, the median BF (Bayes factor; used as a phenotype association score) across the three runs was taken as the output, thereby filtering out spurious results due to stochastic behaviour of the Bayesian computation process. Because preliminary analyses suggested that rare alleles were strongly biasing the results, SNPs with a minor allele count lower than 10 (across all nine populations in the analysed set) were discarded.
As before, bedtools was used to cross-reference the baypass output with our genome annotation to focus on outlier SNPs located in or near genes, yielding a set of genes potentially correlated with CD. As this outlier set was quite large, gene-set enrichment analysis (GSEA) was carried out using topgo version 2.34.0 (Alexa & Rahnenführer, 2016), in order to test for overrepresentation of particular functional gene categories among the outliers. Two GSEAs were performed: one for all genes containing SNPs with a BF score above 15, and one for all genes containing SNPs with a BF score above 20. This increased the chance of identifying larger patterns, and allowed us to gauge the extent to which the GSEA approach was sensitive to the significance cutoff used.

| Testing for sex-linked loci
Because sex-linked variation has been connected to diapause traits in other studies in Lepidoptera (Pruisscher et al., 2020), the combined set of outlier genes from the F ST /π and baypass (BF >20) analyses (n = 379) was tested for Z-linkage. Using blastp (Camacho et al., 2009), the amino acid sequences for the outlier gene models was cross-referenced against a database of all Z-chromosome gene sequences from a recent, chromosome-level assembly of the P. aegeria genome (Lohse et al., 2021 [version 1;awaiting peer review]). blast hits with 70% identity or higher were considered matching genes.

| Scanning for population-specific indels
Because one of the identified outlier genes contained a large, population-specific deletion, a follow-up analysis was conducted in order to assess how common this type of structural variation is in the P. aegeria genome. For this analysis, the read depth-filtered data set was not used, so as not to unnecessarily lose data. Instead, a single mpileup file was generated for all nine populations. Positions where read depth equalled zero for all populations were filtered out, and an indel scan was then run in 100-bp windows using a custom awkbased bash script. Each 100-bp window where a single population had zero reads, while all other populations showed an average read depth of 20 or more across the window, was scored as an indel. Note that using this method, potential indels varying within a population, or shared between more than one population, were not counted.

| Pararge aegeria genome assembly
A total of 148.12 million reads (85% ≥Q30) were generated during sequencing, and using different proportions of this raw read data for assembly greatly affected assembly statistics (Appendix S1: Figure  S2). For the 12 initial linked read assemblies generated, N50 values ranged from 1.8 to 16.9 kbp, while total assembly length ranged from 198.1 to 466.7 Mbp. The 100% and 35% assemblies showed the highest N50 values (16.9 and 14.2 kbp, respectively), as well as the largest number of scaffolds at least 50 kbp in length (n = 778 and 661, respectively). The 35% assembly, however, additionally showed superior gene content and quality, with 74% of busco entries being full length and single copy, and the lowest number of missing orthologues (7.3%). Given these results, the 35% assembly was used as the backbone for subsequent scaffolding, polishing, haplomerging, protein scaffolding and annotation.
The final genome was 479 Mbp in length, and consisted of 26,567 scaffolds (of which 95 scaffolds were >1 Mbp), with a maximum scaffold length of 3,041,493 bp, an N50 of 561 kbp and a total of 10.9% non-ATGC characters. busco analysis of the final genome was much improved, with 88% entries complete and single-copy, 5.7% fragmented and 5.7% missing. A total of 23,567 high-quality gene models were annotated, of which 15,993 had functional annotations, and 33% of the genome was found to contain repetitive content.

| Population history and connectivity analyses
Phylogeographical analysis suggested a deep split between southern versus northern + island P. aegeria populations, with the two Skåne populations as the sister group to all other Swedish populations, and all four northern univoltine populations clustering together (Figure 1b). At least two admixture events could be inferred from the residual allelic variation. First, Gotland appears to derive a large part of its genetic variation (~21%) from a population related to the Stockholm population. Second, Southern Skåne is inferred to have received a smaller amount of genetic variation (~5%) that is shared by all northern and island populations. Subsequent migration edges were not considered geographically likely to reflect true gene flow events, and did not notably improve tree fit, so were not considered further (Appendix S1: Figure S3).

| Regions showing possible selection footprints
All three univoltine/bivoltine comparisons yielded genic regions showing low π, high overall F ST and at least one strongly differentiated

| Regions showing association with critical daylength
Using baypass to test for associations between allele frequency and CD across all nine populations, while controlling for their underlying demographic relationships, yielded 1084 genes containing at least one SNP with a BF score above 15. Of these, 346 genes further contained at least one SNP with a BF score above 20 (Figure 2d). Geneset enrichment analysis revealed few robust patterns among the outlier genes. For genes with a BF above 15, the most significantly enriched GO terms included connections to Toll signalling, limb formation, and responses to fungi and toxins (Table 2). Narrowing the set to genes with a BF above 20 instead yielded a rather different pattern of enriched GO terms, relating to amino acid metabolism and subcellular localization/transport of proteins (Table 3).
One notable result that was reflected in the GSEA at both sensitivity levels was the presence of outlier genes relating to rhythmic processes and circadian function. Five circadian clock genes were found to contain SNPs correlated with CD across the studied popu- timeless2 (BF =16) and cryptochrome (BF =15). This also partly accounts for the emphasis on subcellular protein localization/transport among the enriched GO terms (Table 3), as such transport is part of circadian clock function. Circadian-related GO-term enrichment was further boosted by the presence of a 5-HT (serotonin) receptor among the outliers (BF =17).
One of the outlier regions contained two genes previously described as candidate loci for diapause regulation in P. aegeria: kinesin-associated protein 3 and carnitine o-acyltransferase. Previous work indicates that these two genes are in linkage with timeless (Pruisscher et al., 2018). The highest-scoring SNPs in the baypass output (BF > 50) were found clustered in a 30-kb region containing several genes tied to basic cellular functions: importin (subunit) α, the beta subunit of mitochondrial ATP synthase, the ribosomal protein L24e and the ribosomal export protein NMD3. Overall, the baypass outlier set was dominated by genes showing a cline in allele frequency from north to south, reflecting the correlation between CD and latitude in the studied populations. However, the detailed clinal pattern varied greatly between different outlier genes, with some alleles gently shifting frequencies over latitudes, and others showing a more abrupt sharp north/south cutoff (Figure 3). A list of baypass outlier genes (BF >20) is provided in Appendix S2.
Five genic regions (containing a total of seven genes) appeared as outliers in both the F ST /π and baypass analyses; these showed a high degree of overlap with the candidate regions reported by Pruisscher et al. (2018). Results for these regions are summarized in

TA B L E 2
The 10 most significantly enriched "biological process" GO terms for genes with BF > 15 in the baypass analysis. Total = number of genes assigned a given GO term in the Pararge aegeria genome annotation; Outliers = number of genes assigned the same GO term in the outlier set; p = parent/child-adjusted Fisher p-value total, only 12 outlier loci appeared to be located on the Z chromosome (Appendix S1: Table S4). These included period, which is previously known to be Z-linked in Lepidoptera (Pruisscher et al., 2020), but did not include any of the other loci showing a high F ST /low π signal. Hence, the results of the analyses did not appear to be primarily driven by sex-linked variation.

| Large deletion in timeless
Visual inspection of the alignments revealed an otherwise highcoverage region of ~1500 bp in the circadian gene timeless where no Gotland reads had mapped, indicating a deletion unique to this population ( Figure 4). This deletion spanned exon 15 (out of 16) of TA B L E 3 The 10 most significantly enriched "biological process" GO terms for genes with BF >20 in the baypass analysis. Total = number of genes assigned a given GO term in the Pararge aegeria genome annotation; Outliers = number of genes assigned the same GO term in the outlier set; p = parent/child-adjusted Fisher p-value F I G U R E 3 Patterns of allele frequency for baypass outliers. Each row is one SNP; columns show populations, ordered from north to south. (a) Allele frequencies for genes that received a BF above 20 in the baypass analysis. Only the highest-scored SNP per each gene is shown, and each SNP is only represented once, even if it was located in or near multiple genes. SNPs have been clustered vertically according to similarities across populations; the dendrogram shows the clustering structure. (b) Allele frequencies for all baypass outlier SNPs located in circadian genes, sorted by gene: timeless (tim), period (per), cryptochrome (cry1), timeless2 (tim2) and clockwork orange (cwo)

TA B L E 4
Genic regions that ranked as outliers in both the baypass (BF > 20) and F ST /π analyses. "Comparison" indicates which population contrast returned a given region as having high F ST and low π. "BF" shows the Bayes Factor for the highest-scored SNP in each gene. Genes that were among the 15 candidate loci reported by Pruisscher et al. (2018) are marked with an asterisk (*). Note that the outlier region in contig m232 was large (45 kbp) and contained three genes with high BF scores timeless as well as most of exon 14, shortening the predicted TIM protein sequence by 97 amino acids, or 8% of its length. The deletion was judged unlikely to be a bioinformatic artefact, as reads from all other eight populations mapped to this region (Figure 4). In addition, PCR amplification using a pair of primers that spanned this region yielded a long DNA fragment in two Stockholm individuals, but a short fragment in two Gotland individuals, consistent with a 1.5-kbp deletion in the Gotland samples (Appendix S1: Figure S4).
A genome-wide scan of population-specific indel variation revealed putative deletions at an additional 325 locations in the genome, some of a similar size, but most only one or two hundred bp in size (Appendix S1: Figure S5). In 96 cases, the missing sequence data lay within an exon. The majority of putative deletions were in Sundsvall, although that population may have a somewhat higher risk of false positives, as its read depth was generally lower across the genome. Öland showed deletions at a few loci, and Gotland at one additional locus. However, the deletion in timeless was the largest population-specific indel detected in any population.

| DISCUSS ION
Scandinavian Pararge aegeria constitute a set of relatively recently (i.e., postglacially) immigrated populations showing a certain degree of local adaptation to a climate gradient. Consistent with this, the results presented here suggest a genomic landscape of mild interpopulation differentiation, interspersed with spots of stronger divergence with potential adaptive importance. The outlier searches were conducted with unspecific criteria, the aim being to identify any genic regions potentially associated with selection on voltinism variation. Given this, it is notable that the output of both the selection-footprint scans and phenotype association analyses so prominently featured circadian genes, a group of genes that have previously been implicated in generating variation in life cycle regulation in P. aegeria (Pruisscher et al., 2018) as well as other species (see below). The exonic deletion in timeless was found essentially by luck, as our methods included a visual examination of each detected outlier region, revealing the local drop in read depth (Figure 4). Judging by follow-up scans, such large, population-specific deletions are not particularly common in P. aegeria (Appendix S1: Figure S5). The chance nature of its discovery highlights the risk, inherent to typical SNP-based genomic methods, of overlooking structural variation that may be important to adaptation (Hoban et al., 2016); however, bioinformatic tools do exist that may facilitate the discovery of such variation (Chen et al., 2009;Rausch et al., 2012).

| Circadian gene variation
Of the circadian genes emerging as outliers in the present analyses, timeless, period and cryptochrome are the best characterized. The first two are negative oscillators in the core circadian loop across most insects, including butterflies (Sandrelli et al., 2008;Zhu et al., 2008), while cryptochrome (cry1) codes for the main light receptor that allows the circadian clock to synchronize with the daylight cycle (Yuan et al., 2007;Zhu et al., 2008). Allelic variation in all three of these genes has previously been tied to geographical differences in diapause traits: timeless in pitcher plant mosquitoes (Mathias et al., 2007) and two different Drosophila species Yamada & Yamamoto, 2011); period in corn borer moths (Kozak et al., 2019) and parasitic wasps (Paolucci et al., 2016); and cryptochrome in Drosophila triauraria (Yamada & Yamamoto, 2011).
One of the strongest baypass outliers, containing several variable SNPs (Figure 3b), was clockwork orange (cwo). In Drosophila, the CWO protein acts synergistically with the main period/timeless loop by binding to and repressing several circadian oscillator genes (including tim, per and cwo itself) in the late night, thus forming an essential subloop that helps keep circadian rhytmicity at high amplitude (Kadener et al., 2007;Tomioka & Matsumoto, 2015). While this function has also been inferred in lepidopterans, it has not yet been experimentally confirmed (Brady et al., 2021). Less is known about the function of timeless2 (also known as "mammalian timeless" or timeout), although it appears to help entrain the insect circadian clock to light (Benna et al., 2010), and has been tied to population variation in circadian activity rhythms in rice borer moths (Zhu et al., 2019).
Two other outlier loci merit mention for their connection to this subject. While not a core circadian gene, the baypass outlier casein kinase 1α has been tentatively linked to circadian clock regulation in silk moths (Trang et al., 2006). The highest-scored baypass outlier region (BF ≈ 55) contained the gene importin subunit α. Importins control the transport of proteins into the nucleus; this has been shown to be essential to both fly and mammal clock function for importin β, albeit not for importin α (Lee et al., 2015).
Although a central role for circadian genes in the regulation of diapause in particular and photoperiod-controlled traits in general the circadian clock as such forms a part of the daylength-sensing machinery, as opposed to some circadian genes merely having a pleiotropic role in photoperiodism (Bradshaw & Holzapfel, 2010;Meuti & Denlinger, 2013;Saunders et al., 2004). In either case, the connection between the two time-keeping systems complicates adaptive interpretations for two reasons. First, because different latitudes experience different day/night regimes that may need to be compensated for, circadian genes could potentially vary between populations as a result of direct selection on clock function, with no regard to photoperiodic traits (Bradshaw & Holzapfel, 2010).
Second, clock traits such as circadian period length have themselves been shown to correlate with latitude in various species (Hut et al., 2013;Pivarciova et al., 2016), and in some cases this has been linked to allelic variation in clock genes (Dalla Benetta et al., 2019;Kozak et al., 2019). It is unclear if this type of variation is a side-effect of selection on photoperiodism, or if clock traits themselves are locally adapted to different selective environments (Salmela & Weinig, 2019). In general, the close correlation between latitude and CD across the populations studied here means that, despite controlling for demographic background, the baypass analysis cannot necessarily disentangle the effects of selection on diapause traits from selection on circadian rhythmicity or anything else that may covary with latitude. Functional assays of diapause regulation and/or circadian rhythmicity will be required to understand the adaptive role of the genetic variation found across these populations.
The most drastic novel variation was found in timeless, where Gotland P. aegeria lack 8% of the gene's coding sequence. It is difficult to say what the effect the deletion may have. The lepidopteran circadian clock is relatively well described at the level of whole protein interactions (Brady et al., 2021;Zhu et al., 2008), and experimental data exist on the function of specific parts of a few circadian proteins (Zhang et al., 2017), but to our knowledge the lepidopteran TIM protein has not yet been studied at this level of mechanistic detail. Results from silk moths show that, at least on an amino acid sequence level, both TIM and PER contain many of the functional domains described in Drosophila . These include the cytoplasmic localization domain (CLD) characterized in TIM by Saez and Young (1996), which overlaps partially with the Gotland-specific deletion (Appendix S1: Figure S6). In D. melanogaster, this domain prevents TIM from entering the nucleus; its effect is overridden by dimerization of TIM with PER, which allows the core transcriptional feedback loop of the circadian clock to be closed.
Predicting the functional consequences of circadian gene variation across insect orders is problematic, not least because the insect circadian clock machinery is variable both in which components are present and what roles they play (Yuan et al., 2007). In the monarch (Danaus plexippus), the foremost model species for butterfly circadian function, the core feedback loop is not closed by a TIM/PER dimer, as in Drosophila. Instead, this role is filled by CRY2, a lightinsensitive paralogue of CRY, which forms a trimer with TIM and PER (Zhu et al., 2008). Nonetheless, if TIM-PER interaction sites are conserved, the overlap of the Gotland deletion with the putative silk moth CLD is an indication that this mutation may have some effect on the functioning of the TIM protein. Finally, we note an interesting parallel with a previously known case of a length polymorphism in timeless, whose allele frequency cline also correlates with diapause variation Tauber et al., 2007). However, in that case the protein length variation results from alternative start codons rather than a deletion, and hence affects the N-terminus, the opposite end of the TIM polypeptide from where the Gotland deletion is located.

| Population history and adaptive divergence
Overall genetic differentiation among the studied populations was fairly low, with most comparisons yielding an average genomewide F ST of 0.06-0.08 (Figure 1d). The exception was Öland, which showed relatively strong differentiation from all other populations.
This probably stems from the low allelic diversity of this population ( Figure 1c), which appears to have gone through a demographic bottleneck during its history, as also indicated by its long branch in compatible with, and contribute additional detail to, the picture of genetic diversity in Northern European P. aegeria suggested by earlier results (Tison et al., 2014).
aegeria populations, one of the aims of the present study was to investigate the extent to which selective (local adaptation) versus nonselective processes (genetic history) have contributed to the expressed life cycle patterns. A particularly good indicator of selection is functional differentiation between populations despite a certain degree of gene flow, as it suggests that phenotypic differentiation is continually being maintained by local adaptation (Kawecki & Ebert, 2004). Given that univoltine and bivoltine populations clustered separately in the phylogeographical analyses, with little evidence of recent or ongoing genetic exchange between northern and southern mainland populations, such conclusions cannot as easily be drawn here. However, it should be noted that voltinism is a labile trait on ecological timescales (Altermatt, 2010), generated through a combination of genetic and environmental differences between populations (Lindestad et al., 2019), and has been observed to evolve quickly through modification of photoperiodic plasticity (Bean et al., 2012;Yamanaka et al., 2008). For these reasons, a classical phylogenetic-comparative interpretation with a single transition of voltinism states along the population tree may be misleading.
Assuming a postglacial colonization scenario where a univoltine range margin (as today represented by Sundsvall) tracks the upper limit of the viable climate envelope for the species, it is likely that currently bivoltine P. aegeria populations in this region were univoltine when first established.
Interestingly, the bivoltine Gotland population, which in terms of the candidate SNPs for diapause regulation characterized by Pruisscher et al. (2018) has a distinctly "northern" profile, is here shown to harbour additional, unique variants at two of these same candidate loci (timeless and period). In other words, instead of simply carrying a "northern," high-CD allele of these two genes, this population in fact displays a unique, derived allele of both genes.
If these local variants have a phenotypic effect, this could be interpreted as a result of selection for bivoltinism through a lowered CD following colonization by univoltine genetic stock. It is also worth noting that the Gotland population appears to have an intrinsically longer diapause than neighbouring populations, possibly as an adaptation to partial bivoltinism (Lindestad et al., 2020). Differences in the rate of diapause termination are associated with variation at period in the moth Ostrinia nubilalis (Kozak et al., 2019;Levy et al., 2015), marking the novel circadian gene variation in Gotland as a possible parallel.

| CON CLUS IONS
European populations of Pararge aegeria are a useful model system for studying postglacial migration and climate adaptation (Aalberg Haugen & Gotthard, 2015;Hill et al., 1999;Pateman et al., 2016;Tison et al., 2014). The present results expand on earlier analyses in this species by providing a more fine-scale picture of both neutral and putatively adaptive variation. In a wider perspective, it is especially notable that circadian loci, which appear to play a key role in adaptation to geographical variation in seasonality across a range of insect species, again emerge as potentially under selection.
While any hypothetical effects of the novel mutations described here (either on circadian functioning or on photoperiodic diapause induction) are unknown, these loci stand out as promising targets for further investigation in vivo. Pruisscher for input on bioinformatic methods.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5281/ zenodo.5780591.

DATA AVA I L A B I L I T Y S TAT E M E N T
Bioinformatic scripts have been made available through Zenodo (Lindestad et al., 2021a). Sequence data and genome assembly data have been deposited through EMBL-EBI (Lindestad et al., 2021b).