Contemporary N e estimation using temporally spaced data with linked loci

Abstract The contemporary effective population size Ne is important in many disciplines including population genetics, conservation science and pest management. One of the most popular methods of estimating this quantity uses temporal changes in allele frequency due to genetic drift. A significant assumption of the existing methods is the independence among loci while constructing confidence intervals (CI), which restricts the types of species or genetic data applicable to the methods. Although genetic linkage does not bias point Ne estimates, applying these methods to linked loci can yield unreliable CI that are far too narrow. We extend the current methods to enable the use of many linked loci to produce precise contemporary Ne estimates, while preserving the targeted CI width and coverage. This is achieved by deriving the covariance of changes in allele frequency at linked loci in the face of recombination and sampling errors, such that the extra sampling variance due to between‐locus correlation is properly handled. Extensive simulations are used to verify the new method. We apply the method to two temporally spaced genomic data sets of Anopheles mosquitoes collected from a cluster of villages in Burkina Faso between 2012 and 2014. With over 33,000 linked loci considered, the Ne estimate for Anopheles coluzzii is 9,242 (95% CI 5,702–24,282), and for Anopheles gambiae it is 4,826 (95% CI 3,602–7,353).


| INTRODUC TI ON
Effective population size (N e ) is an important parameter in population genetics. It governs the number of mutants in a population, and hence nucleotide diversity and the number of segregating sites (Charlesworth & Charlesworth, 2010;Wang et al., 2016). It also determines the magnitude of genetic drift, and therefore the stability of allele frequencies over time. Estimates of N e are used in population management, for both conservation and pest control purposes (Lehmann et al., 1998;Waples, 1989). For endangered populations, a certain level of N e needs to be maintained to avoid inbreeding or excessive accumulation of deleterious mutations (Wang et al., 1999). In contrast, if one wishes to reduce the size of a harmful species, N e can act as an indicator to monitor the efficacy of the relevant control measures (Antao et al., 2011). All these applications require robust N e estimates from genetic information.
Depending on the questions of interest, some studies estimate N e for over tens or hundreds of thousands of generations, while some, as in this work, focus on a more contemporary time frame, from the most recent to a few generations ago (Waples, 2005).
Linkage disequilibrium (LD) and temporal changes in allele frequency are two established sources of information to estimate contemporary N e (Luikart et al., 2010). The former utilizes LD signal among unlinked and genetically neutral loci, as genetic drift induces nonzero LD (measured by r 2 ) under finite N e (Hill, 1981). The reason why unlinked loci are used in the LD method is because they provide information about the N e of the parental generation of the samples but little about the population histories further backward in time (Hayes et al., 2003;Waples, 2005). Drift also causes allele frequencies to fluctuate over generations, and hence N e can be estimated through measuring the temporal changes in allele frequency from a collection of primarily neutral and unlinked loci. Within the temporal methods family there are the moment-based F statistics (Jorde & Ryman, 2007;Krimbas & Tsakas, 1971;Nei & Tajima, 1981;Pollak, 1983;Waples, 1989), while more advanced likelihood methods have also been developed (Hui & Burt, 2015;Wang, 2001;Williamson & Slatkin, 1999). These temporal methods estimate the harmonic mean of the N e between the two (or more) sampling events (Waples, 2005). Additionally, the temporal signals from loci under selection (or linked selection) also hold information about N e (Buffalo & Coop, 2019;Khatri, 2016;Wang et al., 2016).
The existing temporal methods for N e estimation were developed with the use of unlinked loci, and hence only provide methods to calculate confidence intervals (CI) for such data. While applying these methods to linked loci could still provide the same N e point estimate, the inference of CI requires the assumption of independence among loci (Hui & Burt, 2015;Wang, 2001). Clearly this assumption does not hold for linked loci, where temporal changes in allele frequency are correlated. Previous attempts to work around the problem include estimating N e along sliding windows or by resampling of loci (Jónás et al., 2016), but neither of these seems to solve the problem directly or make the best use of the available data. With the advance of sequencing technologies, such as restriction site-associated DNA (RADseq) and whole-genome sequencing, obtaining data for large numbers of linked loci is becoming increasingly feasible and affordable. This simultaneously creates the issue of pseudoreplication when computing genomic statistics.
This problem arises when loci in tight linkage or proximity provide correlated information (Patterson et al., 2006;Waples et al., 2020). If not handled properly, it will mislead us to overestimate the amount of information contained in the samples. Extending the temporal methods to appropriately incorporate genetic linkage in CI estimation is therefore essential.

| THEO RY
We will first derive the covariance of the changes in allele frequency for a pair of linked loci under the Wright-Fisher model, and then find the same covariance under the presence of sampling error. Then, we will use this result to approximate the sampling distribution of the temporal F statistic when loci are linked. For tidiness, only key results or equations are displayed in the main text, while the full derivations are given in the Appendix S1 for interested readers.

| Covariance of the changes in allele frequency between two linked loci
Consider a pair of neutral biallelic loci i and j with recombination rate c ij , and let the initial haplotype frequencies be p ij0 . In the next generation, the number of the four haplotypes follow a multinomial distribution with size 2N e , and with probabilities equal to the gametic frequencies after recombination. Let p it and p jt be the allele frequencies of the first allele at the two loci at generation t. The covariance between them is where D ij0 is the LD measure between the two loci in the first temporal sample. Waples (1989), Equation 2 derives the variance of allele frequency due to genetic drift, and Equation 1 above can be viewed as an extension of the same formula for the covariance between two allele frequencies.

| With sampling error
Usually true frequencies can only be estimated through sampling.
Here we consider samples are taken with replacement, or after reproduction. This is equivalent to "sampling plan II" of Waples (1989).
The observed haplotype counts can be modelled by another multinomial distribution. Let x i0 and x j0 be the observed frequencies of the first allele at both loci, and S 0 be the diploid sample size (i.e., 2S 0 haplotypes) at generation 0. The covariance between the two observed frequencies is: Similarly, given the true frequencies and diploid sample size S t at generation t, the covariance between the two observed frequencies x it and x jt is: By applying the total law of covariance, the covariance of x it and x jt given p ij0 .is: The first term of Equation 4 follows from Hill and Robertson (1968). The temporal changes in observed frequencies on the two loci are (x it -x i0 ) and (x jt -x j0 ) respectively. The covariance between them is: Because of linearity this covariance can be further broken down into four other covariances. The two terms cov(x i0 , x jt | p ij0 ) and cov(x it , x j0 | p ij0 ) are 0 under sampling plan II, as samples are taken after reproduction without affecting the contents of the gamete pool (Waples, 1989). Therefore, which is the sum of the two quantities calculated in Equations 2 and 4. These raw changes in observed frequencies need to be normalized across loci. Krimbas and Tsaka (1971) suggest the following standardization: This standardized change in observed frequency δ i , like the raw measure, has mean about 0. Hereafter, we will drop the conditional term p ij0 in our covariances as it becomes cumbersome. Now the covariance between their standardized changes in observed frequency is approximately: where r ij0 is the standardized LD measure between the pair of loci i and j at generation 0 (i.e., the first temporal sample). Note that this expression is remarkably similar to Equation 6, with r ij0 replacing D ij0 for standardization, based on the ratios of expectations. One can use this formula to find the variance of the standardized change in observed frequency at one locus, by substituting c ii = 0 and r ii0 = 1, that is, the covariance with itself: which yields the same equation as in Waples (1989), Equation 7b. We will revisit this equation later as it holds the key to find the point estimate of N e . In addition, we can approximate the correlation between two standardized changes in observed frequency, which is the quotient between Equations 8 and 9:

| Point estimation of F, and its distribution
Consider a more general case with K loci. Let F a be the arithmetic average of K squared standardized changes in observed frequency (Krimbas & Tsaka, 1971): F a is a good estimator for F, the variance of the standardized change in observed frequency (Equation 9). Therefore, its expectation is approximately (Waples, 1989): The point estimate N e can be obtained by solving the above equation with known sample sizes and t. Waples (1989) provides an approximate solution: As mentioned, the expectation (Equation 12) holds regardless of genetic linkage. The main difference is the variance and distribution of F a , and therefore the width of CI, depend heavily on the covariance structure among loci.
Let F a, indep be the F a computed from K independent loci. The classical result suggests KF a, indep ∕F is approximately 2 K distributed, where the subscript in 2 K denotes the degrees of freedom, which is also K (Waples, 1989). The 95% CI for F can be found using the 2.5and 97.5-percentile of 2 K : Similarly, we denote F a, linked as the F a statistic calculated from K linked loci. The distribution of F a, linked should be more dispersed than that from K independent loci. To approximate the distribution of KF a, linked ∕F we let R be a K by K correlation matrix for the standardized changes in observed frequency 1 , 2 , …, K , whose elements are as described in Equation 10. This correlation matrix R, being symmetrical and positively definite, has real positive eigenvalues 1 , 2 , …, K . KF a, linked ∕F is approximately distributed as the sum of K independent random variables: where each Q i is independently and normally distributed with mean 0 and variance i . More details about Q 2 can be found in the Appendix S1.
The closed form of Q 2 is usually not known, but its values can be conveniently computer-generated. The CI for F with linked loci can then be obtained from the empirical quantiles of Q 2 . For instance, For the limiting case of having K independent loci, R is simply an identity matrix, whose eigenvalues are all 1, and thus Q 2 is reduced back to 2 K . While F a is a sum-of-ratios statistic (Jorde & Ryman, 2007), we also introduce its ratio-of-sums counterpart F b : shares the same expectation with F a (and hence provides the same N e ) but with a slightly different variance. To find the CI for F with F b , let W 1∕2 be a diagonal matrix whose (diagonal) elements are √ w i . We then compute the eigenvalues of KW 1∕2 RW 1∕2 and use them to generate Q 2 , the approximate distribution of KF b ∕F. The CI for F can be found as described previously (Equation 16).

| S IMUL ATIONS
To summarize, contemporary N e and its CI can be estimated from linked loci via the following steps: 3. Estimate r 0ij for every pair of loci i and j (see Appendix S1).
4. Calculate the correlation matrix R, using the estimates above, alongside other known parameters t, S 0 , S t and c ij (Equation 10).

Compute the eigenvalues of
Generate the empirical distribution of Q 2 from the eigenvalues Computer simulations were run (see Methods below) to verify the theories behind the two F statistics. The main results are shown in Table 1, with additional results in the Appendix S1. The average F a and F b followed the expectations very closely in all simulation settings, with only a few per cent deviation. However, the bias was exaggerated in the N e scale in some cases, particularly when the sample size to N e ratio was small. The ratio-of-sums F b performed better under such scenarios with much lower N e bias. The standard Recombination frequency was 1e-5 between adjacent bp per generation, and the chromosome length was 1e5 bp. True F (fifth column) is calculated via equation 12 given the parameters. For each combination of parameters the average of F a from 1,000 independent simulation was reported in the sixth column, with the corresponding N e in parentheses (calculated via equation 13). The standard deviations are shown in the seventh column. The next two columns present the proportions of runs (out of 1,000) in which the 95% C.I. cover the true F, using the adjusted Q 2 distributions. Phased and unphased data are treated separately (see SI). The next four columns show the same information for F b , the ratio-of-sum statistic, using the same simulated datasets. For comparison, the coverage of the unadjusted 2 K C.I. (i.e. assuming independence, Equation

14
) is shown in the last column.
deviations of the two estimators were comparable, although those for F a were consistently slightly smaller. Perhaps the more important results are the width and coverage of the 95% CI. Those inferred from the adjusted Q 2 distribution covered the true F in about 95% of the simulations, which accurately reflected the desired confidence level. The adjusted CI worked for both F a and F b statistics, and for both phased and unphased data. For comparison, the last column of Table 1 shows the CI coverage calculated as if loci were independent. These 2 -based CI coverages were all below the targeted level, and in some cases the coverage was as low as 53%. Additional simulations were run to confirm the method's robustness towards different levels of genome-wide recombination (see Methods and Appendix S1), and with more linked loci. In short, the adjusted CI remained accurate with the desired coverage under all the scenarios examined.

| APPLI C ATI ON
We applied the method to two sets of Anopheles gambiae  Note: One set of F a and N e were estimated per chromosome arm (Chrom) initially. For each species two samples with diploid sample sizes S 0 and S t were collected 2 years apart (assumed to correspond to t = 20 generations). Pairwise c ij were calculated from physical distances via Haldane's mapping function, using the published recombination frequency (1.4 centimorgan per megabase; Pombi et al., 2006). The 95% CIs for N e (given in square brackets) were calculated from the adjusted Q 2 , with 50,000 realizations. The last column combines genotypic information from both chromosome arms to provide an overall N e estimate for each species. The same table calculated with F b can be found in the Appendix S1.

F I G U R E 1
The empirical distributions of Q 2 = KF a ∕F for Anopheles coluzzii (left) and Anopheles gambiae (right), where K is the total number of linked loci used on both arms of chromosome 3 (33,148 for A. coluzzii, 33,372 for A. gambiae). The histograms in white are 50,000 realizations of Q 2 based on the eigenvalues of R. The red dotted lines mark the 2.5-and 97.5-percentiles of Q 2 , which are the lower and upper confidence interval for F. The associated upper and lower 95% CI N e are shown in Table 2 (last column). For comparison, the distributions of in the Methods and Appendix S1. With over 33,000 linked loci considered on chromosome arms 3R and 3L combined, the estimated harmonic mean N e in this period via F a were 9,242 (95% CI 5,702-24,282) for A. coluzzii, and 4,826 (95% CI 3,602-7,353) for A. gambiae (Table 2). The empirical distributions of Q 2 for the two estimates are shown in Figure 1. Analyses with F b gave very similar results (Appendix S1).

| DISCUSS ION
Gathering genetic data at two or more time points from a single population ought to be useful for estimating N e , but thus far there have been no methods to estimate proper CI when linked loci are used.
The existing methods implied or assumed independence among loci when inferring CI, as this was the basis for aggregating information across loci (Wang, 2001;Waples, 1989). This assumption severely limits the species applicable to the methods, or the number of loci to be included in one analysis. For instance, existing temporal meth-  (Waples, 1989).
For independent loci, computing var( i ) alone is sufficient to obtain the point and CI estimates (Waples, 1989). For linked loci, however, the now nonzero between-locus covariance cov( i , j ) needs to be considered as drift trajectories are correlated. This work calculates this covariance (Equation 8) from the discrete two-locus twoallele Wright-Fisher model with recombination, and then uses it to provide a method for finding the appropriate CI. The matrix R describes the correlation among i . We can decorrelate them through eigen-decomposition of R. Its eigenvalues are crucial in approximating the sampling distribution of KF a ∕F, and eventually the CI for F and N e . In principle, this eigen-decomposition framework can also help tackle the issue of pseudoreplication on other genomic statistics. One potential example is F ST for population differentiation (Waples et al., 2020). Each locus contributes to the overall F ST but provides duplicated information. Hence var(F ST ) depends on the correlation structure among the individual F ST across loci, which may also depend on parameters such as N e , migration, pairwise LD and recombination rates.
A key component of R is the initial pairwise LD, measured by r ij0 .
During our development we found that using the maximum likelihood estimate (MLE) for r ij0 tends to over-estimate its magnitude and thus the eigenvalues of R, resulting in CIs being slightly more conservative. This is because finite sampling itself also induces some LD (Hill, 1981). Since an unbiased estimator for r ij0 is not found, empirical corrections are imposed, with slightly different treatments for phased and unphased data (see Appendix S1). Another parameter is the pairwise recombination rates c i,j , and a fine-scaled recombination map will be helpful in determining such rates. While a recombination map is not yet available for A. gambiae, we used a published genomewide average recombination frequency (Pombi et al., 2006) and the Haldane mapping function to convert pairwise physical distances into c i,j . Further theoretical work would be needed to investigate the consequences of mis-specifying these rates. While the underlying N e influences the width of CI through affecting the magnitude of average r ij0 (Hill, 1981) and pairwise corr( i , j ), we found that R and the CI are relatively insensitive to mis-specification of N e . A 10-fold over-or under-estimation of N e does not greatly affect the estimated distribution of F (see Appendix S1). Better estimates of R and its components are welcome.
The previous requirement of independent loci precluded the use of high-throughput sequencing technologies, which potentially yield tens of thousands of linked loci. There have been attempts to estimate N e trajectories from whole-genome sequencing data, such as by studying the LD decay curves (Hayes et al., 2003) or identity-bydescent (IBD) tracks (Browning and Browning, 2015), but these onesample methods focus on a much longer horizon (up to thousands of generations) backward in the past. It is also observed that they tend to produce confounded or correlated N e estimates at different time points. Complementary to these works, the temporal method presented here provides an N e estimate that specifically pertains to the time window between samples, without the interference of population size dynamics further in the past (Waples, 2005). It also probes the question of how much genetic information is allocated into the inference of contemporary N e . The LD method suggests that loci with recombination rate c apart contain information about N e of 1/(2c) generations ago (Hayes et al., 2003), which means many loci with shorter distances are not contributing (e.g., c>0.025 is required to infer the most recent 20 generations of N e ). There is a similar claim on the relationship of c and timescale for the IBD method. Our proposed method, in contrast, uses loci with any recombination distances to isolate the drift signal between the sampling events. This temporal method also has the potential to work with RADseq data, where linked loci are discovered from multiple RAD fragments of short length (and hence no IBD information). It is a more accessible alternative to whole-genome sequencing, particularly for nonmodel species, but still generates high-resolution data for demographic inferences (Marandel et al., 2020;O'Loughlin et al., 2014). Note also that the calculation of matrix R only requires LD information from the first temporal sample, which means the second temporal sample could be processed with lower cost technologies such as pooled sequencing, where individual genomes are pooled and sequenced together (Iranmehr et al., 2017;Schlötterer et al., 2014). Thus, costeffective contemporary N e estimation is possible through combining different sequencing technologies with the appropriate experimental designs.
We applied the method to estimate the contemporary N e for A. coluzzii and A. gambiae from a cluster of villages in Burkina Faso.
Their genomes show great diversity, with one variant in about every two bases (Ag1000G, 2020). Although it is possible to use information from multiple chromosomes (see Appendix S1), chromosome 2 was excluded from demographic inferences (Ag1000G, 2020).
Initially, an N e was provided from each of the two chromosome arms 3R and 3L separately. The estimates from both arms were very consistent with overlapping CIs. The genotypes from the two arms were then combined to provide an overall estimate per species, which was more precise. We did not observe significant differences between the estimates from F a and F b because the high number of loci provided sufficient drift signals. However, pseudoreplication was a severe issue here with this amount of loci. For comparison, if one naively treated all loci as independent, then the variance of 2 K would be less than 1/6 of what we calculated from Q 2 (Figure 1). In other words, on average 6+ linked loci provided the same amount of information about genetic drift as one independent locus. Another perspective to the same question is through the "effective number of independent loci" K ′ given a data set (Waples, 2020), which is approximately K � = 2K 2 ∕var(Q 2 ). Note that K ′ is affected by factors such as the level of LD pruning and choice of loci, and hence can only be determined on a case-by-case basis.
Previous N e estimates for Anopheles mosquitos were in the order of 10 6 -10 9 (Ag1000G, 2017; Khatri & Burt, 2019), but those were figures for the entire species rather than a local population, and also averaged over a much longer period backward in time. In contrast, our N e estimates are spatially and temporally restricted. Consistent with previous analyses, we assume the study population is panmictic and closed. If there is immigration into the focal population then this will affect the higher moments of i and hence the N e estimates (Nunney, 2016;Wang & Whitlock, 2003). The window of t = 20 generations is relatively short for immigration to have any significant impact on the local N e estimates (Wang & Whitlock, 2003). A relevant mark-release-recapture experiment shows that the dispersal of Anopheles is mostly short-range (Epopa et al., 2017), and that the migration rate per generation should be small and local. Combining the two factors, our reported estimates should reasonably represent the local and contemporary N e . In principle, it may be possible to jointly infer N e and immigration rates if there are data from more than two time points or multiple geographical locations, but further theoretical work would be needed. Although Nei and Tajima (1981) comment that the effect of selection on temporal F is generally minor, it is a known consequence that selection affects the temporal change in allele frequency (Jónás et al., 2016). For example, mean-reverting balancing selection dampens var( i ), while with directional selection allele frequency changes are positively correlated over time. In our data analyses different measures were introduced to minimize the interference from selection. First, only single nucleotide polymorphisms (SNPs) annotated as intergenic were included. Second, the entire chromosome 2 was excluded because it contains regions with multiple segregating inversions and insecticide-resistance alleles (Ag1000G, 2017; Ag1000G, 2020). Even though there are intergenic SNPs on chromosome 2, they are more likely to be linked with sites under selection. Third, only one SNP was chosen per 1,000-base-pair (bp) window to further reduce the chance of having tight linkages. While the temporal F uses drift signals from neutral loci, there exist other methods which were developed to incorporate selection signals (Buffalo & Coop, 2019;Khatri, 2016). Other standard assumptions for temporal methods (Waples, 1989) apply to our method as well: individuals are random samples from the population; mutation and selection are assumed to be negligible; generations are nonoverlapping; and the number of generations between samples is known without error. In our case we assumed 10 mosquito generations per year, consistent with previous studies (Ag1000G, 2017O'Loughlin et al., 2016). However, other factors, such as the duration of wet and dry seasons, and whether the mosquitoes aestivate during the dry season (Lehmann et al., 2010), may affect this value.
The estimated population size is linearly proportional to t (Equation 13), so that if there were only half as many generations between time points, then the estimated N e would be halved.
Our F a statistic (Equation 11) follows that of Taskas and Krimbas (1971). There are alternative forms of temporal F, such as F c and F k (Waples, 1989), but they are numerically indistinguishable when loci are plentiful with reasonably large sample size and minor allele frequency (MAF) cut-off. Our second statistic F b is a ratio-of-sums statistic similar to that of Jorde and Ryman (2007). In our simulation the mean F a and F b were always within a few per cent from their expectations. This shows that both statistics are good estimators of the rate of loss in genetic diversity due to drift, which is at the scale of 1∕2N e . The bias is unfortunately exaggerated when translated into N e because of the reciprocal relationship, particularly when the denominator term (F a − (1∕2S 0 ) − (1∕2S t )) of Equation 12 is close to 0.
Comparing the two statistics, F b alleviated some bias in cases when the sample size to N e ratio was small, or more generally when there was a lack of drift signals, consistent with previous study (Jorde & Ryman, 2007). F b differs from F a only by the weighing scheme, as F a is the simple arithmetic average while F b is a weighted average of 2 i in which the weights are their heterozygosity. Loci at lower frequencies are given smaller weights and hence contribute less to the overall point estimate. However, there is a trade-off between accuracy and precision, with the unequal weighting system of F b leading to a higher variance (Jorde & Ryman, 2007). Having very few samples and a short sampling horizon may bias F and N e estimates, as investigated previously (Waples, 1989). These are intrinsic problems to the entire temporal F family and unfortunately our methodology is not immune to them. Excluding loci with MAF <5% is a standard safeguard to avoid extreme values in the denominator.
Likelihood-based methods may be better than the moment-based F (Williamson & Slatkin, 1999), and the development of such advanced models for linked loci is worth investigating. Genetic linkage is an important factor affecting the width of CI, on top of sample sizes, K, N e and t (Waples, 1989), and with so many variables affecting the precision we suggest running preliminary simulations to determine appropriate sample sizes. From a practical point of view, we recommend using loci that are sparsely spread along the chromosome with lower average linkage and LD, preferably after some "LD pruning", to lower var(Q 2 ). Although it is tempting to include all available loci for N e estimation, the marginal benefits of adding more linked loci from the same chromosome will diminish, as most information has already been captured. Instead, the excess sampling errors (on allele frequencies, pairwise LD and recombination) and computing burden may outweigh the benefits. Furthermore, it is possible to combine temporal information from more than two time points (Buffalo & Coop, 2019;Williamson & Slatkin, 1999).
The computation effort required by our method is manageable.
Storing R can be memory-hungry when many loci are involved. The calculation of all pairwise r ij0 and c ij0 is repetitive but parallelizable.
Computing the eigenvalues of R is the most computationally demanding task, but a mainstream workstation can easily evaluate R with 33,000+ loci with an optimized linear algebra pack, as demonstrated in our worked example. The maximum number of loci a computer can handle will depend on memory.

| ME THODS
Computer simulations were run using fastsimcoal version 2.6 (Excoffier & Foll, 2011). Two temporally spaced genetic samples separated by t = 10 generations were simulated with chromosome length of 1e5 bp. The mutation rate was 1e-6 per bp per generation, and the recombination frequency was 1e-5 between adjacent bp per generation. Thus, the recombination rate between an arbitrary pair of loci with y bp apart is 0.5 * [1 − (1 − 2 × 10 −5 ) y ]. Different combinations of N e , sample sizes and K were tested. Loci with MAF <5% at either time point were excluded. The two statistics F a and F b were calculated as described. While the simulator outputs phased (haplotypic) data, unphased genotypic data were mimicked by randomly pairing two haplotypes. When calculating initial r ij0 within R, slightly different treatments are required for phased and unphased data (see Appendix S1). In total, 10,000 realizations of Q 2 were generated from the eigenvalues of R (or KW 1∕2 RW 1∕2 for F b ). The 2.5-and 97.5-percentiles were used to calculate the upper and lower CI of F. These calculations were repeated for 1,000 independent simulations for each parameter combination. Results are displayed in Table 1.
Additional simulations were run to confirm the method's robustness towards different levels of genome-wide recombination (as a ratio to chromosome length). Three recombination frequencies (1e-4, 1e-5 and 1e-6 between adjacent bp per generation) with 100-fold differences were chosen to represent cases with high, moderate and low recombination. The chromosome length for all three scenarios were set to 1e5 bp. We also ran another set of simulations to validate the method with K = 20,000 linked loci, a scenario comparable to our real data examples. The results from both additional simulations can be found in the Appendix S1.
The suggested method was also applied to two real data sets of Anopheles mosquitoes (Ag1000G, 2017). Full details about population sampling, sequencing, accessibility, variant calling and filtering can be found in Ag1000G (2020). A brief summary of the procedure is as follows. A total of 290 Anopheles mosquitoes were collected from Burkina Faso at two time points in July 2012 and in July 2014.
July is around the peak of the wet season and of mosquito activity.
Samples were sequenced at the Wellcome Sanger Institute on the Illumina HiSeq200 platform, with targeted coverage of 30 × per individual. Sequence reads were aligned to the AgamP3 reference genome. gatk was used for SNP discovery and filtering. Furthermore, some samples were sequenced a second time using the same method in order to increase the accuracy. The samples went through several stages of quality control to remove samples of poor quality. About 6 million biallelic SNPs on chromosome arms 3R and 3L were annotated as intergenic. Only those with MAF ≥5% at both time points were included. Lastly, to avoid tightly linked SNPs, we randomly chose one SNP per 1,000-bp window, which yielded about 15,000-18,000 SNPs per chromosome arm per species for N e estimation. We provided N e estimates per chromosome arm per species as well as the overall N e estimates per species. N e estimates calculated via F a are reported in Table 2

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest. The funding bodies had no direct role in the design of the study nor in the collection, analysis, interpretation of data or in the writing of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Mosquito data collected in 2012 belong to the Ag1000G Phase 2 release and are publicly available on the MalariaGEN website (https:// www.malar iagen.net/data/ag100 0g-phase -2-ar1). Those collected in 2014 belong to the Ag1000G Phase 3 release (https://www. malar iagen.net/data/ag100 0g-phase 3-snp). The computer codes for data analyses can be found in this public Github repository (https:// github.com/tinyu hui/tempo ral_Ne_linked). The subsets of genotypes used in the study are attached to this paper as a zip file, and are also available in the Github repository above.