Inferring dispersal across a fragmented landscape using reconstructed families in the Glanville fritillary butterfly

Dispersal is important for determining both species ecological processes, such as population viability, and its evolutionary processes, like gene flow and local adaptation. Yet obtaining accurate estimates in the wild through direct observation can be challenging or even impossible, particularly over large spatial and temporal scales. Genotyping many individuals from wild populations can provide detailed inferences about dispersal. We therefore utilized genomewide marker data to estimate dispersal in the classic metapopulation of the Glanville fritillary butterfly (Melitaea cinxia L.), in the Åland Islands in SW Finland. This is an ideal system to test the effectiveness of this approach due to the wealth of information already available covering dispersal across small spatial and temporal scales, but lack of information at larger spatial and temporal scales. We sampled three larvae per larval family group from 3732 groups over a six‐year period and genotyped for 272 SNPs across the genome. We used this empirical data set to reconstruct cases where full‐sibs were detected in different local populations to infer female effective dispersal distance, that is, dispersal events directly contributing to gene flow. On average this was one kilometre, closely matching previous dispersal estimates made using direct observation. To evaluate our power to detect full‐sib families, we performed forward simulations using an individual‐based model constructed and parameterized for the Glanville fritillary metapopulation. Using these simulations, 100% of predicted full‐sibs were correct and over 98% of all true full‐sib pairs were detected. We therefore demonstrate that even in a highly dynamic system with a relatively small number of markers, we can accurately reconstruct full‐sib families and for the first time make inferences on female effective dispersal. This highlights the utility of this approach in systems where it has previously been impossible to obtain accurate estimates of dispersal over both ecological and evolutionary scales.


| INTRODUCTION
The increasing availability of genomic information continues to transform the way we study natural populations. It is now possible to accurately and efficiently measure a wide range of important parameters that directly influence the fitness and survival of wild populations such as effective population size (Gilbert & Whitlock, 2015;Palstra & Fraser, 2012), effective number of breeders (Ackerman et al., 2017), extra pair paternity (Firth, Hadfield, Santure, Slate, & Sheldon, 2015;Griffith, Owens, & Thuman, 2002), heterozygosity Saccheri et al., 1998), inbreeding depression (Huisman, Kruuk, Ellis, Clutton-Brock, & Pemberton, 2016) and reproductive success (Coltman et al., 1999). Another key ecological parameter is dispersal, the ecological and evolutionary causes and consequences of which have been studied for decades. While for some species, including many birds and mammals, very detailed information on individual dispersal events can be obtained by tracking individuals in their natural habitat (e.g., Telfer et al., 2003;Tøttrup et al., 2012), for many others measuring dispersal in the field is difficult or even impossible (Clobert, Danchin, Dhondt, & Nichols, 2001). The use of molecular data makes it possible to examine the role of dispersal and gene flow in influencing the long-term and large-scale spatial genetic structure of populations in such species (Slatkin, 1985), providing an opportunity to investigate dispersal at previously unobtainable ecological and evolutionary scales, as well as in organisms where no dispersal data has previously been available. Additionally, this approach allows the estimation of effective dispersal, that is, dispersal directly associated with reproductive fitness, often a more biologically relevant parameter than general dispersal.
A large number of different estimators have therefore been proposed to infer dispersal from marker data. Initially, F statistics were used to indirectly infer spatial distribution of genetic variation using assumptions from Wright's island model (e.g., Dobzhansky & Wright, 1943). However, it became apparent that this approach has several limitations (reviewed in Whitlock & McCauley, 1999), including the inability to disentangle contemporary versus historical gene flow and dispersal, as well as poor performance in nonequilibrium populations where assumptions of the Wright's island model are violated.
An alternative approach has been to infer dispersal events through first reconstructing relatedness between individuals (e.g., Lepais et al., 2010;Schunter, Pascual, Garza, Raventos, & Macpherson, 2014). For example, if two closely related individuals (e.g., parent-offspring or siblings), whose relationship is established through molecular markers, are sampled in different populations, a dispersal event can be inferred. The utility of this approach has previously been limited by large uncertainties in relatedness estimates resulting from the limited number of available genetic markers and hence difficulties in statistically discriminating even closely related individuals from unrelated. For example, the central idea of relatedness estimation using genetic data is to infer genealogical relationship by comparing similarities in their multilocus genotypes given the Mendelian laws of inheritance. Therefore, the performance of these estimators depends on the number of loci (and alleles at a particular locus), allele frequencies and variance in relatedness within the populations from which the estimates are derived (Csilléry et al., 2006). With the advent of high-throughput sequencing (NGS) technology, it is now both relatively easy and economically feasible to genotype many individuals at hundreds or even thousands of loci, something that will considerably reduce uncertainty surrounding the relatedness estimates. This uncertainty has been further reduced with significant methodological developments, including moving from looking at pairwise comparisons between individuals (e.g., Jones & Ardren, 2003), to considering multiple individuals simultaneously to reconstruct family relationships (Wang, 2007). This not only increases power by taking into account the genotypes of both of the parents, but can also reduce errors when inferring relationships (Wang, 2007), allowing the accurate reconstruction of even pedigrees in wild populations. A major improvement to pedigree reconstruction methods has been to consider the joint likelihood of relatedness among all individuals in the sample (Wang, 2004(Wang, , 2007Wang & Santure, 2009). This method is implemented in the software COLONY (Jones & Wang, 2010), which also considers missing genotypes and genotyping errors (Wang & Santure, 2009), something that can have large influence on the relatedness inferences if not modelled correctly (Wang, 2004). Incorporating error rates is particularly important given that current applications of pedigree reconstruction are likely to use NGS data that can have high error rates and missing genotypes (Glenn, 2011). It is therefore an ideal time to apply these methods to systems where it has previously been challenging or even impossible to have accurate estimates of dispersal across relevant biological spatial scales (Broquet & Petit, 2009).
Here, we report a study on the Glanville fritillary butterfly (Melitaea cinxia), in which molecular markers (SNPs) were used to infer dispersal across a highly dynamic natural metapopulation system over a sixyear period. The Glanville fritillary has been studied for more than two decades in the Åland Islands in Finland, in a large network of about 4,000 small dry meadows (Hanski, 2011). While mark-release-recapture studies of individual butterflies (Hanski, Kuussaari, & Nieminen, 1994;Kuussaari, Nieminen, & Hanski, 1996), tracking studies of individuals with harmonic radar (Ovaskainen et al., 2008) and studies of recolonization of currently unoccupied but suitable habitat patches (van Nouhuys & Hanski, 2002) have provided much information about the extent and spatial scale of dispersal (Hanski, 2011), it has not been possible to examine dispersal at larger temporal and spatial scales. For example, previous work has been unable to directly estimate the levels of gene flow among populations, a parameter of key importance for long-term persistence (Saccheri et al., 1998).
To address this, we sampled around 10,000 individuals (three larvae sampled from each larval nest) from three different regions of the metapopulation over a six-year period and genotyped them for 272 SNPs across the genome. By reconstructing family groups of full-sibs with known geographical position, it is possible to reconstruct the distance an individual female moved between egg laying, thereby inferring breeding dispersal distance. We additionally gained insight into another ecological parameter of importance to the population survival and dynamics of the species, namely the estimated number of surviving larval clutches produced by a single female during her lifetime.
The large geographical and temporal scale of this study provides an unprecedented insight into dispersal movements and performance of individuals in a natural insect system.

| Study species and sampling
The Glanville fritillary butterfly is common in the Åland Islands in SW Finland, where it inhabits a large network of approximately 4,000 dry meadows that have been fully surveyed since 1993 (Hanski, 2011).
Females lay their eggs in large clutches of 150-200 eggs at intervals of one to several days depending on weather conditions (Boggs & Nieminen, 2004). The larvae live gregariously until the last larval instar in the following spring, having spent the winter in a silken "winter nest" spun by the larvae at the base of the host plant. The offspring of two or more females may become mixed either because several females oviposit on the same host plant (M. Saastamoinen, pers. obs.) or because the larval groups merge following short movements while the larvae switch from one host plant individual to another (Kuussaari, van Nouhuys, Hellmann, & Singer, 2004). For further details of the life history, see Ehrlich and Hanski (2004) and for a description of the study area and field methods (Ojanen, Nieminen, Meyke, Pöyry, & Hanski, 2013).
In late August to early September in each year, the entire 4,000 meadows in the patch network are surveyed for the number of larval groups. We estimate that the probability of missing an existing population (false negatives) is around 10%, and these populations are mostly very small, consisting of just one or a few larval groups (Ojanen et al., 2013). The probability of detecting a larval group in a population is estimated to be around 50%-60% (Ojanen et al., 2013).
Since 2007, we have sampled three larvae per larval group for experiments and for DNA sampling. In this study, we analysed data collected from three separate regions in the years from 2007 until 2012 (Fig. S1). Our primary region of interest was the mainland area of Saltvik, where we have sampled a network of 235 patches. For comparison, we also sampled from the islands of Föglö and Sottunga, which have patch networks with 125 and 49 patches, respectively. Sottunga was unoccupied in 1991, at which point 62 larval groups were translocated there from the Finström region of Åland, which neighbours Saltvik, and the metapopulation has persisted ever since.

| DNA extraction and genotyping
Larval tissue was homogenized prior to extraction using TissueLyser (Qiagen) at 30/s for 1.5 mins with Tungsten Carbide Beads, 3 mm (Qiagen). DNA was extracted using the NucleoSpin 96 Tissue Core Kit (Macherey-Nagel). Where DNA yield was low, extracted DNA underwent two rounds of Whole Genome Amplification (WGA) (LGC Genomics). Genotyping was performed using a panel of 272 SNP markers on the KASP platform. Three separate criteria were used to select the panel of SNP markers. Markers from candidate genes relating to flight or dispersal traits, or genes that were differentially expressed after an experimental flight treatment were selected on the basis of previous studies (n = 184; (Saastamoinen, Ikonen, Wong, Lehtonen, & Hanski, 2013;Somervuo et al., 2014;de Jong, Wong, Lehtonen, & Hanski, 2014;Kvist et al., 2015). Putatively neutral SNPs (n = 40) from noncoding regions of the genome were also selected, with the remaining 44 SNPs chosen to ensure that all chromosomes were represented, based on linkage map information (Rastas, Paulin, Hanski, Lehtonen, & Auvinen, 2013). Initial SNP calling was performed on RNA-seq data using the "mpileup" function from "SAMtools" (Li et al., 2009) using default parameter values. This data set included 40 unrelated individuals sampled from across the Åland islands. Only SNPs which had a minor allele frequency (MAF) > 0.2, call rate > 0.9 and SNP quality score > 100 were retained. Neutral SNPs were selected using SOLiD matepair-1 genome sequences (Ahola et al., 2014) and called using an in-house method (Kvist et al., 2015;Rastas et al., 2013). Only heterozygous SNPs spanning all 31 chromosomes from noncoding regions were selected, as the genomic sequences originated from a single male individual. More detailed information on SNP validation and genotyping is reported in . Any SNP or individual with a call rate of <0.95 was removed from further analyses.

| Reconstruction of families
We used COLONY (version 2.0.58 available from https://www. zsl.org/science/software/colony) software (Wang, 2004;Wang & Santure, 2009) to infer relationships (both full-and half-sib) between pairs of sampled individuals using the marker data described above.
COLONY uses a full-pedigree method for sibship inference and parentage assignment (as opposed to pairwise inferences of relatedness) whereby sampled offspring are assigned to hypothetical maternal and paternal families (Wang, 2004;Wang & Santure, 2009). To find a good partition, COLONY uses simulating annealing for optimization, scoring each putative partition using a likelihood-based score (Wang, 2004;Wang & Santure, 2009). Although multiple mating is rare (approximately 6%-8% of the females mate more than once in the wild (Boggs & Nieminen, 2004)), we allowed this in COLONY by assuming polygamous mating in both sexes. We also allowed for inbreeding to occur (Saccheri et al., 1998) and a genotyping error rate of 1.1% per allele (S.C. Wong. pers. comm).
The analyses were carried out separately for each of the six years and the three areas, resulting in 18 different runs. Allele frequencies were estimated on a per year and population basis. Uncertainty in family assignments as a result of the simulated annealing method was reduced by repeating each analysis five times for every data set (i.e., year and population) and retaining only those family relationships that were present in the results from all runs.

| Validating the performance of COLONY using simulated data
The accuracy of sibship reconstruction depends on factors such as sample size, number of sampled individuals for each family and the number and variability of the genetic markers. COLONY assumes that individuals in a sample are taken from a large, randomly mating population (Jones & Wang, 2010). Although our data sets violate this assumption, the accuracy of sibship assignments is not highly sensitive to population genetic structure (Wang & Santure, 2009). To assess the quality of the inferred sibships, we generated simulated data using a forward individual-based model constructed and parameterized for the Glanville fritillary metapopulation in the Åland Islands (see Appendix S1 for the description of the simulation model). We used single-and two-patch models to capture the general population genetic structure in the study areas and test the sensitivity of the COLONY predictions to the level of population genetic structure seen in the Glanville fritillary metapopulation. Although the abundance of larvae widely varied across years, within regions, as well as also across regions (Table S1), for the purpose of this analyses we set the population size near 500 (close to the mean sample size of 462), as computation in COLONY is demanding for large populations. We included all individuals from one generation in both one-patch and two-patch models (492 and 1,014, respectively), as the simulation model effectively randomly sampled at most a handful of offspring from each mother by retaining only offspring that were to survive to adulthood. We finally checked that the simulated data resembled the empirical data in key aspects such as family size distribution, minimum allele frequency, and the number and variation of marker loci. COLONY was then run on these data sets using the same settings as above. The accuracy of inferred sibships from COLONY could then be compared to the simulated data to assess error rates.
COLONY also inferred half-sib families, but these were not consistent across the five replicate runs, with some runs having as many as 48% additional predicted half-sib pairs compared to the consensus results. This result was verified with the simulated data, which demonstrated that 47% of the inferred half-sibs by COLONY were incorrect in the single population data and 44% in the two-patch data set. Therefore, in the remaining analyses, we focused on the full-sib families only.

| Dispersal distances
The yearly spatial scale of dispersal was estimated by fitting dispersal distance distributions to the distances of reconstructed dispersal events from the region of Saltvik (n = 230) assuming either a linear model movement (lm) or a diffusion approximation of a random walk movement model (rw) with constant settling rates (Nathan, Klein, Robledo-Arnuncio, & Revilla, 2012;Turchin, 1998). The parameters were estimated using the probabilistic modelling language Stan (Stan Development Team 2016). The annual parameters determining the scale of dispersal were modelled as lognormally distributed with Half-Cauchy priors of mean 0 and scale 1 on the mean and standard deviation of the lognormal distribution for both models. The probability density functions and the associated mean dispersal distances are given in equations (X1-X4).

Connectivity of individual patches was calculated as
where A is the area of the target patch in hectares, 2a is the mean dispersal distance, d is the distance between patches in kilometres, and N is the number of winter nests in the source patch.

| Spatial genetic structure
As a comparison with the direct estimates of dispersal distances from full-sib reconstruction, we additionally quantified the spatial scale of genetic structure by testing for associations between interindividual relatedness and geographical distance. The rate of decline of interindividual relatedness with distance is expected to reflect the scale of dispersal (Epperson, 2005;Hardy & Vekemans, 1999). We used SPAGeDi 1.5 (Hardy & Vekemans, 2002) to calculate pairwise kinship among individuals (Loiselle, Sork, Nason, & Graham, 1995) in the primary study region of Saltvik for each year. Mean kinship among individuals was calculated within patches, and among individuals within 12 distance classes that were predefined in one-kilometre intervals up to a maximum of 12 km. We used jackknifing over all loci to obtain mean estimates and standard errors of kinship coefficients, and tested the significance of correlations between kinship and geographical distance by permuting spatial locations of individuals 999 times. We randomly subsampled a single individual per larval triplet to attempt to control for family effects. Due to large sample size in 2012, we down sampled the analysis in this year to 120 randomly chosen patches to speed up computations. We repeated this analysis using all patches for all years with only putatively neutral SNPs and found the same pattern.
We additionally calculated F ST for each region and year to assess overall genetic differentiation. Only a single individual per larval triplet was included per patch, and patches with fewer than four families were excluded to prevent bias in estimates due to low sample size (Willing, Dreyer, & van Oosterhout, 2012). Calculations were made using the adegenet package in R (Jombart, 2008; R Core Team 2016).

| Population sizes and structure of empirical data
The number of larval groups in Åland is highly dynamic with large yearly differences in each of the three sampled regions. All three regions saw an initial decline in the number of larval groups, followed by a recovery towards the end of the sampled period (Table S1, Fig. S2).
However, the recovery of Sottunga was far less pronounced than that

| Full-sib family structure across patch networks
On average, over 71% of inferred full-sibs originated from a single larval group, but this varied between populations and years, ranging from 56% in Föglö in 2009 to 90% in Sottunga in 2011 (Table 1).
Out of the larval groups containing at least two individuals (80% of all larval groups, as some sampled individuals either had unsuccessful DNA extractions or were excluded in the SNP filtering process), 76% consisted of full-sibs only, while the remaining larval groups included less related individuals. Table 1 shows the distribution of full-sib families among three categories: those in which all members were located in a single larval group, those in which full-sibs occurred in two or more larval groups but within a single population, and those families in which family members were found in two or more different populations, implying that the female had dispersed between the populations. Temporal differences in the percentage of full-sib families found in different populations were much greater than spatial differences (

| Full-sib family inference consistency and verification using simulations
Predicted full-sib family sizes varied between 1 and 20 (Fig. S3). The five replicate COLONY runs for full-sib inference gave highly consistent results for 17 of the 18 data sets (6 years times 3 regions), with less than 3.3% additional full-sib pairs inferred in a replicate run when compared to the consensus result obtained by accepting only T A B L E 1 The number of larval groups containing two or more individuals (i.e., those included for the subsequent calculations), the predicted number of full-sib families and the proportion of inferred full-sib families on a single plant, within a patch and in multiple patches for each region and year the full-sib pairs inferred in all five runs. In the case where a single population was modelled, 100% of the predicted full-sibs were correct and 98% of all true full-sib pairs were detected by COLONY. In the simulated data based on a two-patch island model with weak genetic differentiation (F ST = 0.03, see Appendix S1), we also found that 100% of the predicted full-sibs were correct and 99% of all true full-sib pairs were detected, demonstrating that COLONY predictions with our marker data are robust to the population genetic structure observed in the study metapopulation.

| Movement patterns
The distance between larval groups belonging to the same full-sib family can be used as a measure of female dispersal distance, as the female must have dispersed between laying the two egg clutches. Figure 1 shows one example of inferred movements between populations, from Saltvik in 2012, while Figure 2 shows the distribution of movement distances. 75% of movement distances are less than 1 km, and 85% are less than 2 km. There are three very long movement distances, more than 6 km. To further validate that these were real dispersal events, we checked the relatedness among the pairs of individuals in the respective families to ascertain that they all supported full-sib relatedness. The estimated scale of mean dispersal distances ranges from 0.66 km to 1.50 km for the negative exponential kernel and 0.74 to 1.48 km for the diffusion kernel (Table 2).

| Spatial genetic structure
We found significant negative relationships between kinship and distance for all years of study, and this was fairly consistent across the years (Figure 3). For 2007, 2008 and 2011, we found significant associations among individuals in distance classes up to two kilometres. In the remaining years, associations were significant up to three kilometres, although mean relatedness among individuals in the three kilometre distance class was quite low (Figure 3). Mean kinship within patches ( Figure 3) and F ST across patches in Saltvik (Table S2) (Table S2).

| DISCUSSION
The ability to reconstruct relatedness of individuals is a powerful tool in the study of ecology and evolution in natural populations.
Here we demonstrate that even in a

| Evaluation of relatedness estimates
As the performance of relatedness estimators depend on many factors, testing them with data generated from simulations with realistic assumptions and a known pedigree can increase the confidence level of the predictions. In addition, the flexibility of simulations allows researchers to test hypotheses that are difficult to do empirically. Despite this advantage, we are aware of only a limited number of studies that analysed empirical sibship data in combination with simulated data (Charman, Sears, Green, & Bourke, 2010;Lepais et al., 2010). Therefore, to assess the reliability of COLONY's predictions for our SNP data sets, we ran COLONY on simulated SNP data sets with known genealogical relationships among individuals. We were primarily interested in knowing whether the number of SNPs (272) was sufficient, given the genetic architecture (i.e., linkage), minimum allele frequency, the distribution of family sizes and the mating system of the species. Although modest numbers of SNP markers (60-120) are reported to be sufficient for accurate reconstruction of sibships (Ackerman et al., 2017;Wang & Santure, 2009), the accuracy can depend on the complexity of genetic structure and other intrinsic population characteristics in a particular data set (Kopps, Kang, Sherwin, & Palsbøll, 2015). Furthermore, our study species violates a few assumptions in COLONY analyses (e.g., no genetic linkage among loci and random mating), although the algorithm is known to be robust to nonrandom mating and linkage of markers (Wang & Santure, 2009).
For the relatively small and well-connected networks of Sottunga and Föglo, a single-patch model was deemed adequate, conforming better to COLONY's assumptions, as moderate gene flow among subpopulations is expected in these networks. For the large network of Saltvik, some spatial genetic structure could be expected especially in years with low abundance (Nair, Fountain, Ikonen, Ojanen, & van Nouhuys, 2016;Orsini, Corander, Alasentie, & Hanski, 2008); hence, we also included a two-patch model with a typical level of genetic differentiation among the subpopulations (semi-independent networks, SINs) (F ST = 0.03, Nair et al., 2016). The results suggest COLONY predictions are not sensitive to the observed level of population genetic differentiation for detecting full-or half-sibs, although the software performed poorly for half-sibs in both the single and two-patch models.
In these analyses, we included all the simulated individuals from one generation for simplicity. We were not concerned with the effects of sampling because sample sizes were relatively large (ranges from 21 to 3,200 larvae) and because three individuals were randomly sampled exhaustively from all the nests. COLONY uses allele frequencies estimated from the data to compute the likelihood, and therefore, more accurate estimates of allele frequencies from larger samples usually T A B L E 2 Observed and estimated annual scale of dispersal in the region of Saltvik F I G U R E 3 Mean kinship and associated standard errors among individuals in Saltvik across distance classes for each year of study. The first distance class (x = −1) shows mean kinship coefficients for within patch comparisons. Filled points above zero indicate individuals were significantly more similar than would be expected due to random sampling of the population lead to more accurate sibship predictions (Wang, 2012). We also did not simulate genotyping errors as they were thought to be fairly low (1.1%; S.C Wong pers. comm., but note that we allowed for genotyping errors in the empirical analyses when inferring full-sibs, see methods).
These simplifications are limitations of our simulated data, but the purpose was to test the power and ability of COLONY to infer sibships using our empirical data set (i.e., with the same number of SNP loci, sample size, population structure).
The results with the simulated data indicated that the predictions from COLONY were highly accurate for full-sib pairs (100% of full-sib families were correctly predicted), while the performance was much poorer for half-sib pairs. This corroborates the results from the empirical SNP data, where predictions were highly consistent for full-sibs but inconsistent for half-sibs among separate COLONY runs. While it was unfortunate not to be able to accurately detect half-sibs, it has been shown that accurately detecting half-sibs is inherently difficult (Ackerman et al., 2017;Wang, 2012;Wang & Santure, 2009).
While it may be tempting to give a recommendation on the minimum number of SNP markers required to accurately estimate fullsib groups, previous work has shown this to be highly dependent on the focal study system. Wang (2012)

| Inferring dispersal parameters from molecular data
In cases where larval groups were found in two different populations, it was possible to infer female effective dispersal distances, with the median dispersal distance for females being 406 m (mean: 1,001 m, range 67-11,004 m). The estimated average dispersal distance is consistent with previous estimates from mark-release-recapture (MRR) data, which range from 280 m to 544 m (Kuussaari et al., 1996;Niitepõld, Mattila, Harrison, & Hanski, 2011). Our results highlight significant annual variation in the scale of dispersal (Table 2), something not possible to detect with traditional MRR approaches, implying a large variation in the scale of colonization and extinction dynamics. Both linear movement and random walk models lead to similar estimates of the scale of yearly movement, which is not surprising given that while biologically implausible, the linear model of movement can lead to similar distribution of dispersal distances compared to the random walk model (Hawkes, 2009;Turchin, 1998). For the first time, we were also able to estimate the frequency and scale of long-distance dispersal events in this system. Our estimate of the longest dispersal event is somewhat higher than those estimated based on MRR studies. Ovaskainen et al. (2008), however, used harmonic radar to track the movement of females butterflies in an area of 30 ha during a two hour period, and the maximum distance moved by an individual was 5,490 m. This shows that while unlikely, the very few long-distance dispersal events estimated in this study, are at least biologically plausible. This is particularly the case as we are potentially observing total dispersal across the entire breeding season with the genetic data, rather than that observed over a two hour period.
Previous estimates of female dispersal rate based on MRR data range from 8% to 40% with great variation among the studies (Hanski et al., 1994;Kuussaari et al., 1996;Niitepõld et al., 2011). Potential reasons explaining the variation in dispersal propensity include variation in weather conditions and habitat configuration among the studies. Factors such as size of a patch, butterfly density, habitat quality and even allele frequency in Pgi genotype have shown to influence dispersal propensity in the Glanville fritillary butterfly (Hanski et al., 1994;Kuussaari et al., 1996;Niitepõld et al., 2011). Similarly, models inferring dispersal from the long-term data on the population dynamics have further estimated that the number of habitat patches visited by a female during her life varies from one to ten (Zheng, Ovaskainen, & Hanski, 2009). Our results on the number of full-sib families found in two or more patches fall on the low side of the estimates on dispersal events, but we are only measuring successful colonization events. Our analysis cannot detect dispersal where the female lays eggs in only a single patch or whether she has moved before she laid her first clutch. However, our data show that female effective dispersal, or colonization distance (estimated here), is similar to the dispersal distance found in MRR data. It is this colonization distance that is the evolutionarily important feature of dispersal, as we are directly capturing gene flow, something not possible with observational studies.
We also used an alternative approach to our dispersal estimates based on family reconstruction, by assessing the spatial decay of kinship coefficients. By contrast to dispersal distances inferred from sibship reconstruction, we detected significant structure up to 3 km.
However, this is expected because spatial genetic structure contains the signal of dispersal and mating that has occurred over many generations. It is also important to note that because kinship within each distance class is compared to mean kinship of the total sample, the measured extent of spatial genetic structure will vary depending on the sample taken. Thus, to truly estimate dispersal from patterns of spatial genetic structure one must control for effective population density, which is particularly difficult to estimate for metapopulations (Hardy et al., 2006;Vekemans & Hardy, 2004). Previous work has also demonstrated a lag time between major demographic events and genetic structure in Åland (Orsini et al., 2008), further highlighting the added value of directly reconstructing dispersal through sibship analysis, which allowed us to track gene flow within single generations. While using family inference to infer the rate and scale of dispersal is not new and has been used to infer dispersal in plants (Dow & Ashley, 1996;Ellstrand & Marshall, 1985), mammals (Peacock & Smith, 1997;Telfer et al., 2003), fish (Jones, Planes, & Thorrold, 2005), birds (Woltmann, Sherry, & Kreiser, 2012) and other insects (Charman et al., 2010;Lepais et al., 2010), in this study we have demonstrated the ability to infer female dispersal even on the large scale of the Åland islands, where the direct observation of so many individuals would be impossible.

| Yearly and regional differences in estimated dispersal
We  (Table 1). This may be explained by the proportion of patch occupancy across years. For example, in both Sottunga and Saltvik, the years preceding the high dispersal years had the lowest proportion of occupied patches recorded across the study period (Table S1, Fig. S2). Following years where patch occupancy is low, the number of colonizations tends to be highest following the reciprocal metapopulation dynamics observed over the lifetime of the annual survey (Hanski, 2011;Ojanen et al., 2013). Negative density-dependent emigration has also been observed in M. cinxia, suggesting that butterflies tend to stay rather than disperse away from high-quality patches, and disperse when local density is low (Kuussaari et al., 1996). This could further explain patterns in Saltvik, where increased dispersal was observed in years following population bottlenecks (i.e., low total number of nests; Table S1, Fig. S2). However, an individual's decision to emigrate will also depend on a number of other factors such as local patch quality and patch boundary effects (Moilanen & Hanski, 1998), and weather during flight season (Kuussaari, Rytteri, Heikkinen, Heliölä, & von Bagh, 2016). Further work is required to determine the relative importance of these factors in explaining regional and yearly variation in dispersal.

| Estimation of other demographic parameters
From the field collected data, the vast majority of collected larval triplets included only full-sibs, while the remainder were most likely mixtures of the offspring of two or more females. This conclusion is drawn based on the following. Firstly, based on previous knowledge only 6%-8% of the females in the wild have been shown to mate more than once (Boggs & Nieminen, 2004). Secondly, of those females that have mated more than once, paternity analyses have indicated that around 80% of the clutches are still sired by only one male, in most cases the first one (Sarhan & Kokko, 2007). Finally, low levels of larval mixing have been reported from previous field studies (Kuussaari et al., 2004). Larval mixing may occur either through two or more females laying eggs on the same host plant or by the larval groups mixing during the summer when they move from the foliated host plant to the next. In the field, approximately 10% of the host plants that contained eggs were observed to have more than one clutch (Kuussaari et al., 2004), but it is unknown whether these clutches were laid by one or more females.
In the majority (71%) of cases, full-sibs occurred in a single larval group. Previous studies, under semi-natural conditions, have shown that some females are able to lay up to ten clutches in their lifetime, with the average number of clutches laid by females being three (Saastamoinen, 2007). In 20% of the cases, the full-sibs occurred in more than two larval groups within a single population and only 9% of the full-sibs occurred in two or more larval groups in more than two local populations, indicating a successful dispersal event (i.e., establishment). Our results here indicate that of those females that successfully reproduce, most are producing only one clutch that survives over the summer. Considering this is crucial when making predictions about a populations survival and persistence.

| CONCLUSIONS
This study provides an example of the power of using relatively small numbers of genetic markers to uncover previously unknown ecological parameters. Caution should still be applied when designing such studies as the number of markers required will be highly dependent on the traits studied, the number of individuals sampled, and the population characteristics (e.g., level of population structure). This study also highlights the continuing value of combining genomic resources with long-term longitudinal field studies, without which many inferences of ecological and evolutionary dynamics in natural populations would be impossible.